Create Equally spaced coordinates on a California State map created with basemap

Question

I have created a California state map with Basemap python library and this shapefile.

Code is below

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(20,20))
map = Basemap(llcrnrlon=-124.48,llcrnrlat=32.53,urcrnrlon=-114.13,urcrnrlat=42.01,
             resolution='c', projection='lcc', lat_0 =  36.778259, lon_0 = -119.417)
#westlimit=-124.48; southlimit=32.53; eastlimit=-114.13; northlimit=42.01
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#f2f2f2',lake_color='#46bcec')
map.drawcoastlines()
map.readshapefile('./CA_Counties/CA_Counties_TIGER', 'CA_Counties_TIGER')
plt.show()

Output

Now what I want to do next is to plot equally spaced points only on the California map which would fill the entire california map (no ocean and nearby states). Example would be the grid shown below as all the points in the grid fills the grid(almost).

The way I think it can be done is by considering California as a polygon or multipolygon and generating equally spaced longitude and latitude inside it. I was looking at https://gis.stackexchange.com/questions/207731/how-to-generate-random-coordinates-in-a-multipolygon-in-python but it didn't really solve my problem as when i ran the following code to generate points

import fiona
from shapely.geometry import shape
import random
from shapely.geometry import Point
from shapely.geometry import Polygon

def generate_random(number, polygon):
    list_of_points = []
    minx, miny, maxx, maxy = polygon.bounds
    counter = 0
    while counter < number:
        pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
        if polygon.contains(pnt):
            list_of_points.append(pnt)
            counter += 1
    return list_of_points
all_points=[]
for pol in fiona.open('./CA_Counties/CA_Counties_TIGER.shp'):
    #l.append(pol['geometry'])
    all_points.append(generate_random(50, Polygon(pol['geometry']['coordinates'])))

It gave me the following error

ValueError: A LinearRing must have at least 3 coordinate tuples

Also I am not even sure if the above code worked would it have been given me all the equally spaced points(lat and lon) which would fill the entire California map. Is there other way to do it or someone can help me with the above code. It must generate lat and lon tho. So i can plot them with map.plot() function in basemap. Also would like to get all the points in a datastructure. For example a list of lists or array of arrays or any other if it works well (maybe a dictionary)

Please identify the specific line which gives you the error, this will help you to learn where the problem is. I don't see how your code could cause the error though, that LinearRing problem looks more like a problem with polygon definition is SHP file. — Michael Entin, Jul 12 '19 at 19:50

score 3 · Accepted Answer · edited Jun 20 '20 at 09:12

I don't fully follow your example, especially why you are using random generation for your point coordinates. If I understand you correctly, you would like to generate points along a regular grid, however only within a given shapefile.

My suggestion is as follows:

import numpy as np
from shapely.geometry import Point
from shapely.ops import cascaded_union
import geopandas as gpd

def generate_grid_in_polygon(spacing, polygon):
    ''' This Function generates evenly spaced points within the given GeoDataFrame.
        The parameter 'spacing' defines the distance between the points in coordinate units. '''
    
    # Convert the GeoDataFrame to a single polygon
    poly_in = cascaded_union([poly for poly in polygon.geometry])

    # Get the bounds of the polygon
    minx, miny, maxx, maxy = poly_in.bounds    
    
    # Now generate the entire grid
    x_coords = list(np.arange(np.floor(minx), int(np.ceil(maxx)), spacing))
    y_coords = list(np.arange(np.floor(miny), int(np.ceil(maxy)), spacing))
    
    grid = [Point(x) for x in zip(np.meshgrid(x_coords, y_coords)[0].flatten(), np.meshgrid(x_coords, y_coords)[1].flatten())]
    
    # Finally only keep the points within the polygon
    list_of_points = [point for point in grid if point.within(poly_in)]

    return list_of_points

A usage example of the function is as follows:

shape_in = gpd.read_file('path/to/shapefile.shp')

points_in_poly = generate_grid_in_polygon(0.25, shape_in)
points_in_poly_gdf = gpd.GeoDataFrame(geometry=points_in_poly)

ax1 = shape_in.plot(facecolor='none', edgecolor='k')
ax1 = points_in_poly_gdf.plot(ax=ax1)

For California with a spacing of 0.25°, this could look like this:

Tweak at your discretion. Depending on the chosen spacing this can take a long time to process. You might want to use spatial indexing for performance improvement.

Also is it possible for you to update the code with spatial indexing. That would be helpful — Kartikeya Sharma, Aug 09 '19 at 18:06

Create Equally spaced coordinates on a California State map created with basemap

1 Answers1