There are many use cases in GIS where we need to compute parameters or aggregate information based on fishnet. One of the use case is to aggregate information for million of points on fishnet. Then use this aggregated fishnet layer to render on map. This aggregated fishnet layer will be more intuitive in nature than million of points.
Source code is given at the end of tutorial.
We will use GeoPandas and Shapely to create fishnet. We will follow following steps to create fishnet
- Import GeoPandas and Shapely
- At line number 7, read shapefile or any other gis dataset using geopandas. In my case, I have one shapefile for which we will create fishnet.
- At line number 10, we are re-projecting the shapefile to projected coordinate system which is 3857. Our original shapefile is in geographic coordinate system that is 4326. The reason to re-project from 4326 to 3857 is to make use of Meters/KM measurements unit to create regular size fishnet.
- At line number 13, we are getting total bounds or layer extent. Total bounds or layer extent are minX-minY and maxX-maxY coordinates that encapsulates the layer within those coordinates.
- At line number 16, we are getting minX, minY, maxX, maxY points. A pictorial, representation for the same is below:
- We will use only minX, minY, maxX, maxY points to create fishnet.
- At line number 23, we have defined square_distance variable. This variable will be use to create regular sized polygons in fishnet based on it’s values. The value in square_distance variable will be in meters.
- From line numbers 20 to 30 is the main logic where we are creating fishnet. We are utilizing only minX, minY, maxX, maxY points and square_distance to create fishnet.
The logic is very simple. We will start from minX-minY point. In every pass, we are making use of square_distance variable to create geometry and then incrementing x,y points (defined at line number 19) by square_distance.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
import geopandas as gpd from shapely import geometry # Create a fishnet if __name__ == '__main__': # Read the shapefile gdf = gpd.read_file('abbottstown_pa_union.shp') # Reproject to projected coordinate system gdf = gdf.to_crs('EPSG:3857') # Get the extent of the shapefile total_bounds = gdf.total_bounds # Get minX, minY, maxX, maxY minX, minY, maxX, maxY = total_bounds # Create a fishnet x, y = (minX, minY) geom_array = [] # Polygon Size square_size = 1000 while y <= maxY: while x <= maxX: geom = geometry.Polygon([(x,y), (x, y+square_size), (x+square_size, y+square_size), (x+square_size, y), (x, y)]) geom_array.append(geom) x += square_size x = minX y += square_size fishnet = gpd.GeoDataFrame(geom_array, columns=['geometry']).set_crs('EPSG:3857') fishnet.to_file('fishnet_grid.shp') |