A 1000 ways to make grids, but here’s mine.
TL;DR — You can create a fishnet grid with exactly equal cells (minus the border cells) by creating a mesh of equally spaced points within a polygon, and then create Thiessen polygons around each mesh point. Complete code is available at the end of this blog.
—
All grids are equal, some are more equal than others.
One of the first steps in analyzing geospatial data is somehow aggregating multiple coordinate data points on a uniform grid over the area you want to do your analysis on. This will help you in calculating different statistics associated with data you are collecting, while also allowing you to further enrich your model by letting you add additional data points to the grid. This is the geospatial equivalent on discretizing a continuous variable or binning data.
If you run a quick Google search (or prompt on chatGTP), you will see many many ways to create a grid, or fishnet grid. This entry is not a thesis on what technique is better or worse, rather an attempt to document a technique I stumbled on while trying to make the right grid for my project. So here it goes.
Selecting and visualizing the grid area.
First, we import all the packages we need for this project.
import numpy as np
import geopandas as gpd
import folium
import warnings
from shapely.prepared import prep
from shapely.geometry import Point,Polygon
from shapely.ops import cascaded_union
from scipy import spatial
warnings.filterwarnings('ignore')Since I live in the city of Karachi in Pakistan, let’s take the examples of a few districts within this city as test case.
# read and filter relevant shapefiles
file_path = '../data/district_shapefiles/District_Boundary.shp'
districts_gpd = gpd.read_file(file_path)
khi_districts_list = ['KARACHI SOUTH','KARACHI EAST','CLIFTON CANTONMENT']
khi_districts_gpd = districts_gpd[districts_gpd.DISTRICT\
.isin(khi_districts_list)]
khi_districts_gpd = khi_districts_gpd.reset_index(drop=True)Next we combine these three districts, Karachi south, east and clifton cantonment into one polygon. This will be total covered area of our grid. We use folium to visualize this
# Combine multiple polygons into one using a cascaded union.
khi_shape = cascaded_union(khi_districts_gpd.geometry)
loc = list(khi_districts_gpd.geometry[0].exterior.coords)[0]
loc = loc[::-1]
map_ = folium.Map(location=loc,zoom_start=12)
map_.add_child(folium.GeoJson(data=khi_shape))
map_
Generate a mesh of points within selected polygon
Since we have the total covered area in a polygon object, the next step is to create a mesh of equidistant points within this polygon.
# Changing projection form WGS84 to Web Mercator for caclulating distance
# between mesh points in meters not radians.
poly = gpd.GeoDataFrame({'geometry':[khi_shape]},crs=4326)\
.to_crs(epsg=3857)['geometry'].values[0]
#Get the bounding points of the selected polygon.
latmin, lonmin, latmax, lonmax = poly.bounds
prep_polygon = prep(poly)
valid_points = []
points = []
# Each point in mesh will be 500 meters apart.
# You can play with this to adjust the cell size
resolution = 500
# create mesh
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))))
# only keep those points that are within the the selected polygon.
valid_points.extend(filter(prep_polygon.contains,points))
# change projection back to WGS84 convert list of points to Geodataframe.
mesh_gpd = gpd.GeoDataFrame({'geometry':valid_points},crs=3857)\
.to_crs(epsg=4326)
mesh_gpd = mesh_gpd.reset_index()Let’s visualize this mesh using folium
map_ = folium.Map(location=loc,zoom_start=12)
for ind, row in mesh_gpd.iterrows():
point = row.geometry
lat, lon = point.y, point.x
folium.CircleMarker([lat, lon], radius=0.1,
color='blue', fill=True,
fill_color='blue').add_to(map_)
map_
Generating Voronoi Tessellations
These points in the mesh will act as the centroids of our Voronoi tessellations, but before we go any further, let’s take a slight detour to explain what are Voronoi Tessellations or Thiessen Polygons.
In mathematical terms, a Voronoi diagram is a partition of a plane into regions close to each of a given set of objects. In the simplest case, these objects are just finitely many points in the plane (called seeds, sites, or generators). For each seed there is a corresponding region, called a Voronoi cell, consisting of all points of the plane closer to that seed than to any other.
Essentially, if a random point spawns inside the region of an existing seed or site, it is the closest the seed of that region. There is an excellent blog on Towards Data Science that dives into this concept in more detail, and I recommend you go read it if you want to learn more. The visualization below has been stolen from the same place.
For the purpose of creating our grid, we need to create Voronoi tessellations around the equidistant mesh of points we created. Because the seeds are equally spaced, the resulting regions or polygons have the same size within the bounds of our selected polygon.
Now let’s run our final bit of code to get our grid
# Getting x and y coordinates from geometry object and using it to
# generate tesselations.
mesh_gpd['lng_lat'] = mesh_gpd.geometry\
.apply(lambda coord: (round(coord.x,5),round(coord.y,5)))
vor = spatial.Voronoi(list(mesh_gpd['lng_lat'].values))
# Function I got off some gist from GitHub that converts infinite
# bordering regions to finite regions.It just works? Apologies to
# the original author, I am unable to locate the original gist and
# give due credit.
def func_voronoi_finite_polygons_2d(vor, radius=None):
"""
Reconstruct infinite voronoi regions in a 2D diagram to finite
regions.
Parameters
----------
vor : Voronoi
Input diagram
radius : float, optional
Distance to 'points at infinity'.
Returns
-------
regions : list of tuples
Indices of vertices in each revised Voronoi regions.
vertices : list of tuples
Coordinates for revised Voronoi vertices. Same as coordinates
of input vertices, with 'points at infinity' appended to the
end.
"""
if vor.points.shape[1] != 2:
raise ValueError("Requires 2D input")
new_regions = []
new_vertices = vor.vertices.tolist()
center = vor.points.mean(axis=0)
if radius is None:
radius = vor.points.ptp().max()*2
# Construct a map containing all ridges for a given point
all_ridges = {}
for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
all_ridges.setdefault(p1, []).append((p2, v1, v2))
all_ridges.setdefault(p2, []).append((p1, v1, v2))
# Reconstruct infinite regions
for p1, region in enumerate(vor.point_region):
vertices = vor.regions[region]
if all(v >= 0 for v in vertices):
# finite region
new_regions.append(vertices)
continue
# reconstruct a non-finite region
ridges = all_ridges[p1]
new_region = [v for v in vertices if v >= 0]
for p2, v1, v2 in ridges:
if v2 < 0:
v1, v2 = v2, v1
if v1 >= 0:
# finite ridge: already in the region
continue
# Compute the missing endpoint of an infinite ridge
t = vor.points[p2] - vor.points[p1] # tangent
t /= np.linalg.norm(t)
n = np.array([-t[1], t[0]]) # normal
midpoint = vor.points[[p1, p2]].mean(axis=0)
direction = np.sign(np.dot(midpoint - center, n)) * n
far_point = vor.vertices[v2] + direction * radius
new_region.append(len(new_vertices))
new_vertices.append(far_point.tolist())
# sort region counterclockwise
vs = np.asarray([new_vertices[v] for v in new_region])
c = vs.mean(axis=0)
angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
new_region = np.array(new_region)[np.argsort(angles)]
# finish
new_regions.append(new_region.tolist())
return new_regions, np.asarray(new_vertices)
# Convering infinite border regions to finite regions.
regions,vertices=func_voronoi_finite_polygons_2d(vor,radius=0.07)
# Funciton that clips our resulting polygon within the
# bounds of our original polygon.
def func_intersect_vor(row):
return khi_shape.intersection(row.vor)
# Converting regions and vertices to Polygon objects
# and getting a geodataframe with the final grid.
polys = []
for region in regions:
poly = vertices[region]
polys.append(Polygon(poly))
vor_gpd = gpd.GeoDataFrame(polys)
vor_gpd = vor_gpd.rename(columns={
0:'vor'
})
vor_gpd['bounded_vor'] = vor_gpd.apply(func_intersect_vor,1)
vor_gpd = gpd.GeoDataFrame(vor_gpd,geometry='bounded_vor')
vor_gpd = vor_gpd.drop(columns='vor')
vor_gpd = vor_gpd.reset_index()
vor_gpd = vor_gpd.rename(columns={'index':'vor_id'})
vor_gpd['vor_id'] = vor_gpd.vor_id.apply(lambda x: 'vor'+str(x))Our grid is ready. Let’s visualize it using folium.
map_ = folium.Map(location=loc,zoom_start=12)
for ind, row in vor_gpd.iterrows():
geometry_obj = row.bounded_vor
map_.add_child(folium.GeoJson(geometry_obj))
map_The vor_gpd variable here is a GeoDataFrame where each polygon object in the geometry column is a cell of the fishnet grid, or one region of the corresponding seed of the Voroni Tesselation — depending on how you want to look it.
All the code and data required to replicate this experiment is in the github repository below. Do you think there’s a better way to do this? (I’m sure there is) Please let me know in the comments section.
