pyshp generate random points in shapefile

Question

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'm having trouble installing `geopandas` in my environment, but yeah that seems like the least-worst option now — Ralph Asher, Mar 04 '20 at 20:41
I suggest using `anaconda` that manages the installation of `geopandas` easily. Let me check with `pyshp`... — dzang, Mar 04 '20 at 20:44

score 0 · Accepted Answer · answered Mar 04 '20 at 21:14

0

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')

answered Mar 04 '20 at 21:14

dzang

2,160
2
12
21

This is awesome! How to convert `gs` to a pandas df for export? – Ralph Asher Mar 04 '20 at 21:26
You can save it as a shapefile with `gs.to_file('points.shp')` or how would like to export it? If you found the answer good, please upvote it and mark it as 'accepted answer' by clicking on gray tick mark on the left of the upvotes number, thanks. – dzang Mar 04 '20 at 21:29
I'd like to save the random points as a CSV with x/y columns so I assume converting to a pandas dataframe first is the right step – Ralph Asher Mar 04 '20 at 21:37
you can modify the for loop changing `points.append({'x': x, 'y': y})` and then create a dataframe `df = pd.DataFrame(points)` which can be saved as a csv in the format you want. – dzang Mar 05 '20 at 08:36

pyshp generate random points in shapefile

1 Answers1