我有一个很大的栅格文件,想要根据多边形或shapefile将其裁剪成6个numpy数组。我有shapefile和多边形作为geopandas数据帧。请问有谁能帮助我使用Python(不使用arcpy)完成此操作?
我有一个很大的栅格文件,想要根据多边形或shapefile将其裁剪成6个numpy数组。我有shapefile和多边形作为geopandas数据帧。请问有谁能帮助我使用Python(不使用arcpy)完成此操作?
import gdal
import ogr, osr
# converts coordinates to index
def bbox2ix(bbox,gt):
xo = int(round((bbox[0] - gt[0])/gt[1]))
yo = int(round((gt[3] - bbox[3])/gt[1]))
xd = int(round((bbox[1] - bbox[0])/gt[1]))
yd = int(round((bbox[3] - bbox[2])/gt[1]))
return(xo,yo,xd,yd)
def rasclip(ras,shp):
ds = gdal.Open(ras)
gt = ds.GetGeoTransform()
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shp, 0)
layer = dataSource.GetLayer()
for feature in layer:
xo,yo,xd,yd = bbox2ix(feature.GetGeometryRef().GetEnvelope(),gt)
arr = ds.ReadAsArray(xo,yo,xd,yd)
yield arr
layer.ResetReading()
ds = None
dataSource = None
shapefile.shp
,栅格图像为big_raster.tif
,您可以按以下方式使用它:gen = rasclip('big_raster.tif','shapefile.shp')
# manually with next
clip = next(gen)
## some processing or inspection here
# clip with next feature
clip = next(gen)
# or with iteration
for clip in gen:
## apply stuff to clip
pass # remove
arr = ds.ReadAsArray(xo,yo,xd,yd)
将读取要素和栅格之间的交集数组。您可以在“yield”之前对“arr”进行任何修改。 - Val