# Creating polygon grid from point grid using Geopandas

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.

``````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.

• As per the help center please do not include chit chat like statements of appreciation within your posts.
– PolyGeo
Oct 23, 2021 at 6:05

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

• 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
– M.O.
Oct 24, 2021 at 2:56

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

``````

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, lon`pair 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
``````
• Please explain what your finding was and what you changed to make it work. Oct 23, 2021 at 13:50