0

When given a shapefile as below as an Input, I would like to generate a 100x100 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

0

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

| improve this answer | |
  • I made use of the code by I end up getting a single polygon , I have attached the images in the questions above. – roshualine Nov 29 '19 at 13:40
  • @roshualine added edit, you are using geographic coordinate system, so it is in degrees (width and height) change it to equal units – Jan Doležal Nov 29 '19 at 15:10
  • Thank you for pointing out my mistake. I was unaware of the conversion part, the link was helpful. – roshualine Dec 2 '19 at 8:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.