There are many use cases in GIS world, where the information has to be aggregated, an easy way to achieve this is via gridding or binning, where the area of interest is divided into small sections called grids or bins.

These sections are mostly of rectangular form (which can be easily converted into geotiffs), but in some cases even circles or hexagons are also used. You can read a good tutorial from mapbox using Qgis with a mmqgis plugin here.

While MMQGIS is a good enough for most use cases, I was looking for a way to create fishnets or grids by giving the bbox of the region and width/height in meters, which would then make it easy to script the process for any dynamic use.

Luckily there is a “recipe” available which takes input in EPSG:3857 projections and hence takes the width/height of grid cell in meter, here.

The code is pretty easy to comprehend and creates simple polygons for each computed cell in the layer.

``````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)
poly = ogr.Geometry(ogr.wkbPolygon)

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

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

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

# Close DataSources
outDataSource.Destroy()

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, sys.argv, sys.argv, sys.argv, sys.argv, sys.argv, sys.argv )
``````

Here is an example usage: python grid.py grid.shp 992325.66 1484723.41 494849.32 781786.14 10000 10000

Here is an example of a 10km/10 km grid Python scripts are generally usable in most use cases, how ever incase one needs a Postgis alternative, i.e a function which runs in Postgis database and creates a fishnet in a table, then It might be useful to use a query which I found here.

``````CREATE OR REPLACE FUNCTION public.makegrid\_2d (
bound\_polygon public.geometry,
grid\_step integer,
metric\_srid integer = 3857 --metric SRID
)
RETURNS public.geometry AS
\$body\$
DECLARE
BoundM public.geometry; --Bound polygon transformed to the metric projection (with metric\_srid SRID)
Xmin DOUBLE PRECISION;
Xmax DOUBLE PRECISION;
Ymax DOUBLE PRECISION;
X DOUBLE PRECISION;
Y DOUBLE PRECISION;
sectors public.geometry[];
i INTEGER;
BEGIN
BoundM := ST\_Transform(\$1, \$3); --From WGS84 (SRID 4326) to the metric projection, to operate with step in meters
Xmin := ST\_XMin(BoundM);
Xmax := ST\_XMax(BoundM);
Ymax := ST\_YMax(BoundM);

Y := ST\_YMin(BoundM); --current sector's corner coordinate
i := -1;
<> LOOP IF (Y > Ymax) THEN --Better if generating polygons exceeds the bound for one step. You always can crop the result. But if not you may get not quite correct data for outbound polygons (e.g. if you calculate frequency per sector) EXIT; END IF; X := Xmin; <> LOOP IF (X > Xmax) THEN EXIT; END IF; i := i + 1; sectors[i] := ST\_GeomFromText('POLYGON(('||X||' '||Y||', '||(X+\$2)||' '||Y||', '||(X+\$2)||' '||(Y+\$2)||', '||X||' '||(Y+\$2)||', '||X||' '||Y||'))', \$3); X := X + \$2; END LOOP xloop; Y := Y + \$2; END LOOP yloop; RETURN ST\_Transform(ST\_Collect(sectors), ST\_SRID(\$1)); END; \$body\$ LANGUAGE 'plpgsql';
``````

I hope this helps.