I want to generate 1000 random points within a specific ZIP Code Tabulation Area shapefile using pyshp. My code is:
import shapefile
zctashape = shapefile.Reader('C:/mypath/tl_2019_us_zcta510.shp')
shapefile_len = len(zctashape.shapes())
#identify the index for zcta 84049 and designate number of points.
zcta_to_use = '84049'
pointcount = 1000
for i in range(0,shapefile_len):
if zctashape.record(i)[1] == zcta_to_use:
record_i=i
break
record_i
How do I proceed from here? I apologize if this is basic, I am mostly an R user and it is very easy in that language. Shapefile downloaded from https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2019&layergroup=ZIP+Code+Tabulation+Areas
I am using geopandas
since I think it makes it easier. I don't know if there is an easy way to achieve this solely with pyshp
:
import numpy as np
import random
from shapely.geometry import Point
import geopandas as gpd
import matplotlib.pyplot as plt
filename = 'tl_2019_us_zcta510.shp'
gdf = gpd.read_file(filename)
#identify the index for zcta 84049 and designate number of points.
zcta_to_use = '84049'
pointcount = 100
aoi = gdf[gdf['ZCTA5CE10'] == zcta_to_use]
aoi_geom = aoi.unary_union
# find area bounds
bounds = aoi_geom.bounds
xmin, ymin, xmax, ymax = bounds
xext = xmax - xmin
yext = ymax - ymin
points = []
while len(points) < pointcount:
# generate a random x and y
x = xmin + random.random() * xext
y = ymin + random.random() * yext
p = Point(x, y)
if aoi_geom.contains(p): # check if point is inside geometry
points.append(p)
gs = gpd.GeoSeries(points)
fig, ax = plt.subplots()
aoi.plot(ax=ax, facecolor='none', edgecolor='steelblue')
gs.plot(ax=ax, color='r')