Convert Latitude/Longitude points to Grid Polygons in GeoPandas

Question

I'm trying to plot data from a text file (organized with latitude, longitude, and pollen flux values) as a raster grid in Python. I'm using the code for Choropleth Map on https://autogis-site.readthedocs.io/en/latest/notebooks/L5/02_interactive-map-folium.html to try to display the data. My GeoPandas geodataframe has point geometries; however, it looks like the geometry of the points in the tutorial are already multipolygons, which I assume are the squares in the grid. How do I convert my data (assuming each latitude/longitude point is the center of a pixel in a grid) into gridded geopandas (geodataframe) data? The projection I'll be using is Lambert Conformal Conic projection.

To clarify what my geodataframe looks like, when doing gdf.head(10).to_dict(), it looks like this

    {'geoid': {0: '0',
      1: '1',
      2: '2',
      3: '3',
      4: '4',
      5: '5',
      6: '6',
      7: '7',
      8: '8',
      9: '9'},
     'geometry': {0: <shapely.geometry.point.Point at 0x7fa3e7feee90>,
      1: <shapely.geometry.point.Point at 0x7fa3e7feed10>,
      2: <shapely.geometry.point.Point at 0x7fa3e7feef90>,
      3: <shapely.geometry.point.Point at 0x7fa3e7fe4f90>,
      4: <shapely.geometry.point.Point at 0x7fa3e7fe4e50>,
      5: <shapely.geometry.point.Point at 0x7fa3e7fe4bd0>,
      6: <shapely.geometry.point.Point at 0x7fa3e7fe4ed0>,
      7: <shapely.geometry.point.Point at 0x7fa3e7fe4c90>,
      8: <shapely.geometry.point.Point at 0x7fa3e7fe4d50>,
      9: <shapely.geometry.point.Point at 0x7fa3e7fe4c10>},
     'pollenflux': {0: 0.0,
      1: 0.0,
      2: 0.0,
      3: 0.0,
      4: 0.0,
      5: 0.0,
      6: 0.0,
      7: 0.0,
      8: 0.0,
      9: 0.0}}

when it should be formatted like this:

    {'geoid': {0: '0',
      1: '1',
      2: '2',
      3: '3',
      4: '4',
      5: '5',
      6: '6',
      7: '7',
      8: '8',
      9: '9'},
     'geometry': {0: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363f50>,
      1: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363c90>,
      2: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e93631d0>,
      3: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363f10>,
      4: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363410>,
      5: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363a90>,
      6: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363d90>,
      7: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363d10>,
      8: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363390>,
      9: <shapely.geometry.multipolygon.MultiPolygon at 0x7fa3e9363190>},
     'pop18': {0: 108,
      1: 273,
      2: 239,
      3: 202,
      4: 261,
      5: 236,
      6: 121,
      7: 196,
      8: 397,
      9: 230}}

Answer 1

I believe the issue you're having is that the code expects a GeoDataFrame with Polygons or MultiPolygons but yours only has Points.

Here's a quick way to generate a new GeoDataFrame with squares around your points:

import shapely
import numpy 

def get_square_around_point(point_geom, delta_size=0.0005):
    
    point_coords = np.array(point_geom.coords[0])

    c1 = point_coords + [-delta_size,-delta_size]
    c2 = point_coords + [-delta_size,+delta_size]
    c3 = point_coords + [+delta_size,+delta_size]
    c4 = point_coords + [+delta_size,-delta_size]
    
    square_geom = shapely.geometry.Polygon([c1,c2,c3,c4])
    
    return square_geom

def get_gdf_with_squares(gdf_with_points, delta_size=0.0005):
    gdf_squares = gdf_with_points.copy()
    gdf_squares['geometry'] = (gdf_with_points['geometry']
                               .apply(get_square_around_point, 
                                      delta_size))
    
    return gdf_squares

# This last command actually executes the two functions above. 
gdf_squares = get_gdf_with_squares(gdf, delta_size=0.0005)

Note that the delta_size parameter dictates the distance between the coordinates of the corners of the squares and the center. When your original data is in WGS84 coordinates (EPSG 4326), using a delta of 0.0005 will result in squares of about 100 meters if your data is in Central Texas.

Take a look at your input data, find the CRS it's using and try to estimate a good delta value that will generate big enough squares but won't overlap with each other.

Hopefully that gets the rest of the code working.