Python Geopandas – Creating Polygon Grid from Point Grid

geopandaspythonvector-grid

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.

Best Answer

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

Related Question