根据shapefile或多边形将栅格剪切成NumPy数组

Question

我有一个很大的栅格文件，想要根据多边形或shapefile将其裁剪成6个numpy数组。我有shapefile和多边形作为geopandas数据帧。请问有谁能帮助我使用Python（不使用arcpy）完成此操作？

- Lisa Mathew

---

你看过py-gdal cookbook吗？ - Val

是的，我的 shapefile 中有 6 个不同的多边形，并且我想将我的栅格文件分成 6 个不同的 numpy 数组，每个数组对应 shapefile 中的一个多边形。我该怎么做呢？ - Lisa Mathew

假设你的多边形是独立的要素，你可以遍历要素，然后依次使用每个要素将文件剪裁。 - Val

谢谢。这些多边形是独立的要素。我不知道如何迭代要素以剪切大栅格文件。您能否详细说明如何剪切？ - Lisa Mathew

Answer 1

我创建了一个小型生成器，可以满足你的需求。我选择使用生成器而不是直接迭代功能，因为如果您想检查数组，这更方便。如果你愿意，你仍然可以迭代生成器。

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

- Val

---

很高兴能够帮到您。如果解决方案回答了您的问题，请考虑接受它。 - Val

我做了。再次感谢。最后一个问题。在创建单独的数组时，是否可以排除无数据像素？ - Lisa Mathew

如何排除？行arr = ds.ReadAsArray(xo,yo,xd,yd)将读取要素和栅格之间的交集数组。您可以在“yield”之前对“arr”进行任何修改。 - Val

如何获取图像中所有像素值的坐标？我正在使用以下链接作为参考：https://gis.stackexchange.com/questions/260304/extract-raster-values-within-shapefile-with-pygeoprocessing-or-gdal。 - Lisa Mathew

假设有一个矩形多边形的shapefile，如果它们是不规则的呢？ - Ricardo Guerreiro