creating fishet grid using python

Question

I'm trying to create a fishnet grid to aggregate some geospatial data. I am using the following code from the GDAL/OGR Python cookbook however it returns only one polygon object which is basically a huge rectangle. Please let me know what I'm doing wrong.

source: https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-fishnet-grid

import os, sys
import ogr
from math import ceil


def main(outputGridfn,xmin,xmax,ymin,ymax,gridHeight,gridWidth):

    # convert sys.argv to float
    xmin = float(xmin)
    xmax = float(xmax)
    ymin = float(ymin)
    ymax = float(ymax)
    gridWidth = float(gridWidth)
    gridHeight = float(gridHeight)

    # get rows
    rows = ceil((ymax-ymin)/gridHeight)
    # get columns
    cols = ceil((xmax-xmin)/gridWidth)

    # start grid cell envelope
    ringXleftOrigin = xmin
    ringXrightOrigin = xmin + gridWidth
    ringYtopOrigin = ymax
    ringYbottomOrigin = ymax-gridHeight

    # create output file
    outDriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputGridfn):
        os.remove(outputGridfn)
    outDataSource = outDriver.CreateDataSource(outputGridfn)
    outLayer = outDataSource.CreateLayer(outputGridfn,geom_type=ogr.wkbPolygon )
    featureDefn = outLayer.GetLayerDefn()

    # create grid cells
    countcols = 0
    while countcols < cols:
        countcols += 1

        # reset envelope for rows
        ringYtop = ringYtopOrigin
        ringYbottom =ringYbottomOrigin
        countrows = 0

        while countrows < rows:
            countrows += 1
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            # add new geom to layer
            outFeature = ogr.Feature(featureDefn)
            outFeature.SetGeometry(poly)
            outLayer.CreateFeature(outFeature)
            outFeature = None

            # new envelope for next poly
            ringYtop = ringYtop - gridHeight
            ringYbottom = ringYbottom - gridHeight

        # new envelope for next poly
        ringXleftOrigin = ringXleftOrigin + gridWidth
        ringXrightOrigin = ringXrightOrigin + gridWidth

    # Save and close DataSources
    outDataSource = None


if __name__ == "__main__":

    #
    # example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
    #

    if len( sys.argv ) != 8:
        print "[ ERROR ] you must supply seven arguments: output-shapefile-name.shp xmin xmax ymin ymax gridHeight gridWidth"
        sys.exit( 1 )

    main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7] )

Output:

Answer 1

Found a rather derpy way to do it using voronoi tessellations. Function below accepts a shapely.geometry.Polygon object, use the resolution parameter to control the size of grid cell in meters.

Have tested for multiple locations in Pakistan; Karachi, Lahore, Gujranwala.

from functools import partial
from shapely import ops, geometry, wkt
from shapely.prepared import prep
from pyproj import Proj,transform
import pyproj 
import geopandas as gpd


def get_grid_from_poly(poly,resolution,base_proj='epsg:4326'):
    """
    Takes a polygon and creates a mesh of equally spaced points within it.
    ---
    poly: bounding box polygon within which we want to create mesh of Points
    
    resolution: Distance in meters between each point of the equally spaced mesh of Points.
    This determines the size of each grid cell in the fishnet grid. generate_fishnet function
    esentially creates voronoi tesselations with each point of the mesh treated as a centroid.
    
    base_proj: Base Projection from which coordinates are projected to Web Mercator, EPSG: 3857
    """

    project = partial(
        pyproj.transform,
        pyproj.Proj(init=base_proj),
        pyproj.Proj(init='epsg:3857')
    )

    poly = ops.transform(project,poly)

    # shp = districts_gpd[districts_gpd.DISTRICT==district_name]
    # shp = shp.to_crs(epsg=3857)
    # poly = shp.geometry.values[0]

    latmin, lonmin, latmax, lonmax = poly.bounds

    # construct a rectangular mesh

    prep_polygon = prep(poly)
    valid_points = []
    points = []

    #0.012 degrees is approximately equal to 1km
    # resolution = 1000

    for lat in np.arange(latmin, latmax, resolution):
        for lon in np.arange(lonmin, lonmax, resolution):
            points.append(geometry.Point((round(lat,4), round(lon,4))))

    # validate if each point falls inside shape using
    # the prepared polygon
    valid_points.extend(filter(prep_polygon.contains, points))


    mesh =  gpd.GeoDataFrame(valid_points)
    mesh = mesh.reset_index()
    mesh = mesh.rename(columns={
        'index':'p_id',
        0:'geometry'
        })

    mesh['p_id'] = mesh.p_id.apply(lambda x: 'p{id}'.format(id=str(x)))
    #print(mesh)
    mesh['geometry'] = mesh['geometry'].apply(lambda x: ops.transform(func_flip,x))
    mesh['geometry'] = mesh['geometry'].apply(lambda x: ops.transform(func_flip,x))
    
    mesh.crs = 'EPSG:3857'
    mesh = mesh.to_crs(crs=4326)
    mesh['coords']=mesh['geometry'].apply(lambda x: list(x.coords)[0])
    mesh['lon'],mesh['lat']=mesh.coords.apply(lambda x: x[0]),mesh.coords.apply(lambda x: x[1])
    mesh = mesh.drop(columns=['coords'])
    mesh['lon_lat'] =  mesh.apply(func_lonlat_from_lat_lon,1)
    mesh = gpd.GeoDataFrame(mesh,geometry='geometry',crs='EPSG:4326')
    
    return mesh

Answer 2

The code look good to me, perhaps it's a "flushing" issue. Try also closing the layer by setting it to None:

outLayer = None

The code you are using is quite elaborate, a nice demonstration of how OGR works, but you could simply it quite a bit. Consider the example below:

Define the grid and output:

# settings
output_shp = r'D:\tmp_grid_v5.shp'

drv = ogr.GetDriverByName('ESRI Shapefile')

if os.path.exists(output_shp):
    # DeleteDataSource will delete related files, compared to os.remove
    drv.DeleteDataSource(output_shp)

# grid definition
# extent
ulx, uly, lrx, lry = -180, 90, 180, -90

# resolution 1 degree grid
xres = 1
yres = -1

# half the resolution
dx = xres/2
dy = yres/2

Calculate the center coordinates of each polygon. Using Numpy meshgrid is convenient but calculates all at once. If memory is an issue you could do this one by one.

# center coordinates
xx, yy = np.meshgrid(
    np.arange(ulx+dx, lrx+dx, xres), 
    np.arange(uly+dy, lry+dy, yres),
)

Initialize the output shapefile:

ds = drv.CreateDataSource(output_shp)
lyr = ds.CreateLayer(output_shp, geom_type=ogr.wkbPolygon)
fdefn = lyr.GetLayerDefn()

Loop over each center coordinate and add the polygon to the output:

for x,y in zip(xx.ravel(), yy.ravel()):

    poly_wkt = f'POLYGON (({x-dx} {y-dy}, {x+dx} {y-dy}, {x+dx} {y+dy}, {x-dx} {y+dy}, {x-dx} {y-dy}))'

    ft = ogr.Feature(fdefn)
    ft.SetGeometry(ogr.CreateGeometryFromWkt(poly_wkt))
    lyr.CreateFeature(ft)
    ft = None

Close the output

lyr = None
ds = None