1

我正在尝试获取一个包含 5 个多边形的 shapefile，将它们中的每一个栅格化为一个栅格文件，然后将所有这些栅格相加。我一直在尝试在 python 中使用rasterioand来完成此操作geocube，但遇到了一些问题，即我的光栅文件未对齐，从而阻止我执行求和。我将在下面详细说明我的所有数据和步骤。

数据源：我正在使用这个包含纽约市 5 个行政区的 shapefile：https ://data.cityofnewyork.us/City-Government/Borough-Boundaries/tqmj-j8zm

首先，我想为每个多边形/自治市分配一个虚构的示例“test_value”，这可以是任何东西。我想首先栅格化每个多边形，以便为每个隔离的栅格化多边形创建一个单独的栅格文件。这意味着我将有一个仅对曼哈顿进行栅格化的栅格文件，就好像它是城市中唯一的行政区一样，然后将下一个栅格文件仅对布鲁克林进行栅格化，就好像它是城市中唯一的行政区等。这是我用来完成此操作的代码：

polygons = gpd.read_file('Boroughs_Test/Boroughs.shp')
polygon_IDs = polygons['ID'].tolist()

for i in polygon_IDs:
    x = polygons.loc[polygons['ID'] == i]
    vector_fn = x
    out_grid = make_geocube(
        vector_data=vector_fn,
        measurements=["test_value"],
        resolution=(-25, 25),
        fill=-9999,
    )
    out_grid["test_value"].rio.to_raster(str(i) + "_Output_Raster.tif")

这行得通，我会收到每个行政区的单独光栅文件。例如这里是曼哈顿： 在此处输入图像描述

这是布朗克斯： 在此处输入图像描述

然后我一次查看所有 5 个栅格以构成整个城市： 在此处输入图像描述

因此，现在每个栅格化多边形都有自己独特的“test_value”值。我想在所有栅格中添加所有这些“test_value”值。所以这就是我试图用来总结所有这些栅格层的方法：

for raster in folder:
    blank_raster + raster
    with rasterio.open('Result_raster_NYC.tif', 'w', **result_profile) as dst:
        dst.write(result_array, indexes=[1])

本质上，我遍历目录中的每个自治市镇栅格，并不断将它们添加到全为 0 的空白栅格中。但是，我怀疑这将是正确的方法，并且我相信我遗漏了一些需要用于光栅计算的重要考虑因素。

如何扩充我的代码，以便对目录中的所有栅格图层求和？如果需要指定，这些栅格图层只有一个波段包含“test_values”。最终产品看起来与上面的完整纽约地图相同，但每个栅格化的自治市镇都包含其相应“test_value”值的像素值，具有不同的灰色阴影。

4

0 回答 0