2

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.

1
  • 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

2 Answers 2

3

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

1
  • 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
1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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