python - 根据 shapefile 或多边形将栅格剪辑成 numpy 数组

我有一个很大的光栅文件，我想根据多边形或 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 被称为 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

关于python - 根据 shapefile 或多边形将栅格剪辑成 numpy 数组，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/50035578/