Fastest way to produce a grid of points that fall within a polygon or shape?

I am using shapely in python and trying to generate evenly spaced points in a grid that fall within a shape in the fastest O(n) time. The shape may be any closed polygon, not just a square or circle. My current approach is:

1. Find min/max y & x to build a rectangle.
2. Build a grid of points given a spacing parameter (resolution)
3. Verify one-by-one if the points fall within the shape.

Is there a faster way to do this?

# determine maximum edges
polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# construct a rectangular mesh
points = []
for lat in np.arange(latmin, latmax, resolution):
    for lon in np.arange(lonmin, lonmax, resolution):
        points.append(Point((round(lat,4), round(lon,4))))

# validate if each point falls inside shape
valid_points.extend([i for i in points if polygon.contains(i)])

I saw that you answered your question (and seems to be happy with using intersection) but also note that shapely (and the underlying geos library) have prepared geometries for more efficient batch operations on some predicates (contains, contains_properly, covers, and intersects). See Prepared geometry operations.

Adapted from the code in your question, it could be used like so:

from shapely.prepared import prep

# determine maximum edges
polygon = shape(geojson['features'][i]['geometry'])
latmin, lonmin, latmax, lonmax = polygon.bounds

# create prepared polygon
prep_polygon = prep(polygon)

# construct a rectangular mesh
points = []
for lat in np.arange(latmin, latmax, resolution):
    for lon in np.arange(lonmin, lonmax, resolution):
        points.append(Point((round(lat,4), round(lon,4))))

# validate if each point falls inside shape using
# the prepared polygon
valid_points.extend(filter(prep_polygon.contains, points))