Python – Generating Grid Polygons by Converting Meter to Degree

geopandaspythonpython 3shapefilevector-grid

When given a shapefile as below as an Input, I would like to generate a 100×100 meter polygon grid using python.

I went through Creating a regular polygon grid over a spatial extent, rotated by a given angle and Creating grid polygons from coordinates using R or Python but still am not able to achieve the end result.

Based on Creating polygon grid using Geopandas my code is as below (similar to the post), grids are not created within polygon extent but is created somewhere outside.

Which part have I made a mistake?

import geopandas as gpd
from shapely.geometry import Polygon
import numpy as np
poly = gpd.read_file('point.shp')
xmin,ymin,xmax,ymax = poly.total_bounds
print(xmin,ymin,xmax,ymax)
#-118.77194444444446 38.61038888888889 -118.66961111111112 38.78205555555555
length = 100
wide = 100
cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), length))
rows.reverse()
polygons = []
for x in cols:
    for y in rows:
        polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y-length), (x, y-length)]) )

grid = gpd.GeoDataFrame({'geometry':polygons})
grid.to_file('grid.shp')

shpfile
result Grid

Best Answer

I suggest you to use Geopandas, it is self explanatory in the answer below:

https://gis.stackexchange.com/a/316460/111567

Edit: You are using WGS84 with Lat Lon, so width and height must be in degrees. Right now 100 x 100 means 100 degrees, that is why it is one polygon

Change that to roughly 0.001 (110 km * 0.001 is circa 100 m) for lat (length/Y direction) and wide to 0.00125 (80 km * 0.00125 is circa 100 m) for lon (width/X direction)

import geopandas as gpd
from shapely.geometry import Polygon
import numpy as np
poly = gpd.read_file('point.shp')
xmin,ymin,xmax,ymax = poly.total_bounds
print(xmin,ymin,xmax,ymax)
#-118.77194444444446 38.61038888888889 -118.66961111111112 38.78205555555555
length = 0.001 ##### here
wide = 0.00125 ##### here
cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), length))
rows.reverse()
polygons = []
for x in cols:
    for y in rows:
        polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y-length), (x, y-length)]) )

grid = gpd.GeoDataFrame({'geometry':polygons})
grid.to_file('grid.shp')

Why those numbers? Check for more information about conversion of Lat/Lon to meters: https://www.thoughtco.com/degree-of-latitude-and-longitude-distance-4070616

Related Question