Creating polygon grid from point grid using Geopandas

Question

I have a grid of points, that look something like this, the points represent the bottom left corner of an area and I want to retrieve that grid. I've tried doing this and it's not working, the resulting grid has overlapping squares when it shouldn't have.

original point grid and resulting polygon grid that is wrong

def transform_to_grid(mesh, origin="bottom-left"):
 mesh.sort_values(by=["i_lat", "i_lon"], inplace=True)
 xmax=mesh["i_lon"].max()
 ymax=mesh["i_lat"].max()
 mesh_max=mesh[(mesh["i_lat"]==ymax) & (mesh["i_lon"]==xmax)]
 mesh_min=mesh[(mesh["i_lat"]==0) & (mesh["i_lon"]==0)]
 xmin=0
 ymin=0
 ## i_lat increases from south to north.
 if mesh_max["lat"].values[0]>mesh_min["lat"].values[0]:
 Bool_StoN=True
 ysteps=np.arange(ymin, ymax, 1)
 
 else:
 Bool_StoN=False
 ysteps=np.arange(ymax, ymin, -1)
 ## i_lon increases from left to right.
 if mesh_max["lon"].values[0]>mesh_min["lon"].values[0]:
 BoolLtoR=True
 xsteps=np.arange(xmin, xmax, 1)
 else:
 BoolLtoR=False
 xsteps=np.arange(xmax, xmin, -1)
 
 center=[]
 grid_cells = []
 info0 = mesh[(mesh["i_lat"]==ysteps[0]) & (mesh["i_lon"]==xsteps[0])]
 x0 = info0.lon.values[0]
 y0 = info0.lat.values[0]
 
 for x in range(1, xmax):
 for y in range(1, ymax):
 ## Bounds
 info=mesh[(mesh["i_lat"]==ysteps[y]) & (mesh["i_lon"]==xsteps[x])]
 x1 = info.lon.values[0]
 y1 = info.lat.values[0]
 ## Append information
 grid_cells.append( box(x0, y0, x1, y1) )
 if origin=="bottom-left":
 center.append((ysteps[y-1], xsteps[x-1]))
 elif origin=="bottom-right":
 center.append((ysteps[y-1], xsteps[x]))
 elif origin=="top-right":
 center.append((ysteps[y], xsteps[x]))
 elif origin=="top-left":
 center.append((ysteps[y], xsteps[x-1]))
 y0=y1
 x0=x1
 cell = gpd.GeoDataFrame({"geometry":grid_cells, "indexes":center}, crs="epsg:4326")
 return cell.merge(mesh[["lat", "lon", "indexes"]], on="indexes", how="inner")
grid=transform_to_grid(mesh, origin="bottom-left")
mesh=
 lat lon i_lat i_lon geometry 
0 -31.765354 -71.461792 0 0 POINT (-71.46179 -31.76535) 
1 -31.767292 -71.149109 0 1 POINT (-71.14911 -31.76729) 
2 -31.768250 -70.836365 0 2 POINT (-70.83636 -31.76825) 
3 -31.768250 -70.523621 0 3 POINT (-70.52362 -31.76825) 
4 -31.767292 -70.210907 0 4 POINT (-70.21091 -31.76729) 
.. ... ... ... ... ... 
79 -30.167404 -69.302917 6 7 POINT (-69.30292 -30.16740) 
80 -30.162682 -68.996948 6 8 POINT (-68.99695 -30.16268) 
81 -30.157040 -68.691040 6 9 POINT (-68.69104 -30.15704) 
82 -30.150452 -68.385193 6 10 POINT (-68.38519 -30.15045) 
83 -30.142902 -68.079407 6 11 POINT (-68.07941 -30.14290) 

I am very confused as to why the code doesn't work. I am going from south to north and west to east in the loop, keeping the older lat,lon pair to create the polygon.

What I need to do is basically connect the points (that are not equally spaced) and retrieve the square grid from there.

• PolyGeo: As per the help center please do not include chit chat like statements of appreciation within your posts.

Answer 1

This is a fairly efficient implementation, using (fully vectorized) numpy operations:

# %% Creation, just to have a runnable example
import numpy as np
import geopandas as gpd
import pygeos
# Setup points
x = np.arange(11.0)
y = np.arange(11.0)
yy, xx = np.meshgrid(y, x, indexing="ij") 
points = pygeos.creation.points(xx.ravel(), yy.ravel())
gdf = gpd.GeoDataFrame(geometry=points)
gdf.plot()
# %% The actual work:
# Grab the coordinates
# Reshape into the mesh order
# Build the polygons
# coords = pygeos.get_coordinates(gdf.geometry)
# This might not work, alternatively, convert from shapely first:
coords = pygeos.get_coordinates(pygeos.from_shapely(gdf.geometry))
nrow = 10
ncol = 10
n = nrow * ncol
nvertex = (nrow + 1) * (ncol + 1) 
assert len(coords) == nvertex
# Make sure the coordinates are ordered into grid form:
x = coords[:, 0]
y = coords[:, 1]
order = np.lexsort((x, y))
x = x[order].reshape((nrow + 1, ncol + 1))
y = y[order].reshape((nrow + 1, ncol + 1))
# Setup the indexers
left = lower = slice(None, -1)
upper = right = slice(1, None)
corners = [
 [lower, left],
 [lower, right],
 [upper, right],
 [upper, left],
]
# Allocate output array
xy = np.empty((n, 4, 2))
# Set the vertices
for i, (rows, cols) in enumerate(corners):
 xy[:, i, 0] = x[rows, cols].ravel()
 xy[:, i, 1] = y[rows, cols].ravel()
# Create geodataframe and plot result
mesh_geometry = pygeos.creation.polygons(xy)
mesh_gdf = gpd.GeoDataFrame(geometry=mesh_geometry)
mesh_gdf.plot(color="none")

I'm using a lexsort the get the proper ordering (by rows and columns): this might go wrong if your points are located on a very curvilinear grid. In that case, you'd probably need to project to a Cartesian coordinate system first (with e.g. pyproj).

• M.O.: Thanks! Your version is 28 times more efficient than mine! I was able to get rid of a lot of code that I didn't need based on what was in the data frame. For example, the i_lat and i_lon variables are what tells me what row and column they are so the managing of the data is more straightforward. Hope you see the edit to my answer and let me know what you think and if it can be improved even further

Answer 2

I was able to make it work. For anyone that needs the same (because no matter how much I googled I couldn't find an example online) here is the modified function. It goes S to N and W to E by indexes, it retrieves all four corners by using the value of index in the loop and the one before it. With that it starts filling up the lists to then create the dataframe.

######################################################
def transform_to_grid(mesh, origin=""):
 mesh.sort_values(by=["i_lat", "i_lon"], inplace=True)
 xmax=mesh["i_lon"].max()
 ymax=mesh["i_lat"].max()
 mesh_max=mesh[(mesh["i_lat"]==ymax) & (mesh["i_lon"]==xmax)]
 mesh_min=mesh[(mesh["i_lat"]==0) & (mesh["i_lon"]==0)]
 xmin=0
 ymin=0
 ## i_lat increases from south to north.
 if mesh_max["lat"].values[0]>mesh_min["lat"].values[0]:
 Bool_StoN=True
 ysteps=np.arange(ymin, ymax+1, 1)
 
 else:
 Bool_StoN=False
 ysteps=np.arange(ymax+1, ymin, -1)
 ## i_lon increases from left to right.
 if mesh_max["lon"].values[0]>mesh_min["lon"].values[0]:
 BoolLtoR=True
 xsteps=np.arange(xmin, xmax+1, 1)
 else:
 BoolLtoR=False
 xsteps=np.arange(xmax+1, xmin, -1)
 
 center=[]
 grid_cells = []
 x0=0
 ## Creates polygons from S to N, and W to E
 for x1 in range(1, xmax+1):
 y0=0
 for y1 in range(1, ymax+1):
 ## Bounds
 infobl=mesh[(mesh["i_lat"]==ysteps[y0]) & (mesh["i_lon"]==xsteps[x0])]
 infobr=mesh[(mesh["i_lat"]==ysteps[y0]) & (mesh["i_lon"]==xsteps[x1])]
 infotr=mesh[(mesh["i_lat"]==ysteps[y1]) & (mesh["i_lon"]==xsteps[x1])]
 infotl=mesh[(mesh["i_lat"]==ysteps[y1]) & (mesh["i_lon"]==xsteps[x0])]
 lontr = infotr.lon.values[0]
 lattr = infotr.lat.values[0]
 lonbl = infobl.lon.values[0]
 latbl = infobl.lat.values[0]
 lonbr = infobr.lon.values[0]
 latbr = infobr.lat.values[0]
 lontl = infotl.lon.values[0]
 lattl = infotl.lat.values[0]
 ## Append information
 grid_cells.append( Polygon([[lonbl, latbl], [lonbr, latbr], [lontr, lattr], [lontl, lattl]]))
 if origin=="bottom-left":
 lats.append(latbl)
 lons.append(lonbl)
 indexes.append((ysteps[y0], xsteps[x0]))
 elif origin=="bottom-right":
 lats.append(latbr)
 lons.append(lonbr)
 indexes.append((ysteps[y0], xsteps[x1]))
 elif origin=="top-right":
 lats.append(lattr)
 lons.append(lontr)
 indexes.append((ysteps[y1], xsteps[x1]))
 elif origin=="top-left":
 lats.append(lattl)
 lons.append(lontl)
 indexes.append((ysteps[y1], xsteps[x0]))
 y0=y1
 
 x0=x1
 new_mesh = gpd.GeoDataFrame({"geometry":grid_cells, "indexes":indexes, "Cent_lat":lats, "Cent_lon":lons}, crs="epsg:4326")
 return new_mesh

enter image description here

I am sure it can be done more efficiently, so I am open to suggestions. The errors were that the y0=0 was missing before the loop, the use of Shapely box instead of Polygon, keeping the lat, lonpair of the previous point instead of their indexes, and not retrieving all four corners.

---

Edit:

Based on Huite Bootsma I was able to modify their code to fit my data and it's 28 times more efficient than what I came up with!! That's amazing. Here is the edited version of the code they gave me.

def transform_to_grid_vec(mesh, origin=""):
 mesh.sort_values(by=["i_lat", "i_lon"], inplace=True)
 xmax=mesh["i_lon"].max()
 ymax=mesh["i_lat"].max()
 #######
 y=np.reshape(mesh["lat"].values, (-1, xmax+1))
 x=np.reshape(mesh["lon"].values, (-1, xmax+1))
 # Setup the indexers
 left = lower = slice(None, -1)
 upper = right = slice(1, None)
 corners = [ [lower, left], [lower, right], [upper, right], [upper, left] ]
 n=xmax*ymax
 # Allocate output array
 xy = np.empty((n, 4, 2))
 # Set the vertices
 for i, (rows, cols) in enumerate(corners):
 xy[:, i, 0] = x[rows, cols].ravel()
 xy[:, i, 1] = y[rows, cols].ravel()
 if origin=="bottom-left":
 origin=0
 origin_x=0
 origin_y=0
 elif origin=="bottom-right":
 origin=1
 origin_x=1
 origin_y=0
 elif origin=="top-right":
 origin=2
 origin_x=1
 origin_y=1
 elif origin=="top-left":
 origin=3
 origin_x=0
 origin_y=1
 else:
 origin=0
 origin_x=0
 origin_y=0
 # Create geodataframe and plot result
 mesh_geometry = pygeos.creation.polygons(xy) 
 mesh_gdf = gpd.GeoDataFrame({"i_lon":np.tile(np.arange(origin_x, origin_x+xmax), ymax), "i_lat":np.repeat(np.arange(origin_y, origin_y+ymax), xmax)}, geometry=mesh_geometry, crs="epsg:4326")
 mesh_gdf[["Cent_lon", "Cent_lat"]] = gpd.GeoDataFrame(geometry=mesh_gdf.boundary).apply(lambda x: x.iloc[0].coords[origin], axis=1, result_type="expand")
 mesh_gdf["Area_model"]=mesh_gdf.area
 return mesh_gdf

• bugmenot123: Please explain what your finding was and what you changed to make it work.