2

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 enter image description here

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).enter image description here

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 Generating random coordinates in 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)

11
  • Hi, could you elaborate a bit on the meaning of "randomly generated and equally spaced points" aspect? Do you mean that the placement of the first point should be random or that for a given location it would random if a point is created or not? or that the randomness lies in an attribute of the point (since I see points of different colors) ?
    – vagvaf
    Jul 12, 2019 at 5:38
  • In general, if the points are equally spaced, they are not randomly placed.
    – vagvaf
    Jul 12, 2019 at 5:46
  • Yes I agree the meaning of randomness is confusing here and i think the thing I am trying to do is to create equally spaced points on the California map.So i can plot these points(lat and lon) later on. You can take the randomness away from the scope of this question. Also I would like to get all the points in a data structure like list, dictionary, array etc. It doesn't really matter a lot. Jul 12, 2019 at 5:49
  • 4
    Feel free to edit your question in order to get the right message across
    – vagvaf
    Jul 12, 2019 at 5:51
  • @vagvaf I have edited my question. Thank You! Jul 12, 2019 at 5:56

2 Answers 2

3

I tried this code and it worked

import geopandas as gp
from shapely.geometry import Point

metersteps=10
planningarea=gp.read_file('california.shp')

planningarea=planningarea.to_crs(epsg=32633) # you can leave this out or replace with any other
envelope=planningarea.geometry.envelope
envelopecoords=list(envelope.geometry.exterior[0].coords)

urtuple=envelopecoords[2]
lltuple=envelopecoords[0]


xur=urtuple[0]
yur=urtuple[1]
xll=lltuple[0]
yll=lltuple[1]

regularpoints = gp.GeoDataFrame(columns=['geometry'])

for i in range(round((xur-xll)/metersteps)):
    yll2=yll

    for i in range(round((yur-yll)/metersteps)):
        regularpoint=gp.GeoDataFrame([Point(xur,yll2)],columns=['geometry'])
        regularpoints=regularpoints.append(regularpoint)
        yll2=yll2+metersteps
        print(len(regularpoints))

    xur = xur - metersteps

regularpoints.to_file('regularpoints.shp', driver='Shapefile')

It could be that the order of the upper left and the lower right are not the same in your case. You will have to do a print out of the envelope list and check them yourself.

If you want the points only to be in California just clip the result.

3
  • Thanks for than answer! It gives me all the points but how can I plot these points on the California map using basemap as these are not in lat/lon coordinates. Do you think is it possible to convert these points to lat/lon or obtain these points in lat/lon values. Thank You! Jul 14, 2019 at 1:14
  • yes sure you can convert them. just use this command: regularpoints=regularpoints.to_crs(epsg=4236)
    – Leo
    Jul 15, 2019 at 7:28
  • If you like my answer, mark it with the green arrow as being right :)
    – Leo
    Jul 15, 2019 at 7:28
0

Not sure if you've found a solution for this or not but the verified answer from this question pretty much sums it up. You create a grid of equally spaced points based on min, max for lat and lon and then keep only those points that are inside the polygon in question. If you want to look at the state of California without county lines you could look at doing a cascaded union to get one polygon for the state.

In order to work with this you'll need to make sure you're coordinates are projected to EPSG:3857 or WGS 84 / Pseudo-Mercator

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)

#distance in meters between each grid point/coordinate
resolution=1000

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

Have added code here for better visibility. Hope this helps.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.