1

I have tried to create a fishnet(10 by 10) using the below code. But the output is as shown below. What is that I'm missing in the below code?

import geopandas as gpd
from shapely.geometry import Polygon
import numpy as np
ras= rasterio.open(r"My Raster file")

xmin,ymin,xmax,ymax = ras.bounds

length = 1000
wide = 1000

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")
df= gpd.read_file('grid.shp')
df.plot()

enter image description here

2
  • can you add what the bounds of your raster are? – Ian Turton Jan 21 at 10:04
  • @IanTurton Here is my bounding box BoundingBox(left=80.85505119895836, bottom=16.31986397086808, right=81.85798771751823, top=17.489998966915074) – hashalluring Jan 21 at 10:06
2

You are mixing degrees and metres and it will end badly.

cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), length))

will "expand" to

cols = list(range(int(np.floor(80.85505)), int(np.ceil(81.8579)), 1000))
rows = list(range(int(np.floor(16.31986)), int(np.ceil(17.4899)), 1000))

which is clearly not going to work, you will end up with one box 1000 degrees to the side (as you see in your map).

If you want to work in metres then you will need to switch to a local projected CRS.

2

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.