计算多个重叠多边形（shapefile）内的光栅景观比例（百分比）？

Question

我认为最简单的方法是提取每个多边形内的栅格值并计算比例。是否可以在不将整个网格读取为数组的情况下这样做？

---

我有 1992 年至 2015 年的 23 个年度全球分类栅格（分辨率 = 0.00277778 度）和一个具有 354 个形状的多边形矢量（在某些部分重叠）。由于重叠（自相交），将它们作为栅格使用并不容易。两者都投影在“+proj=longlat +datum=WGS84 +no_defs”中。

栅格由 10 - 220 类组成 该多边形的 ABC_ID 范围为 1 - 449

一年时间看起来像：

---

我需要创建一个表，例如：

---

我已经尝试通过以下方式实现这一目标：

1. 区域统计
2. Pk 工具（从光栅中提取矢量样本）
3. LecoS（叠加栅格指标）
4. SAGA GIS 的“交叉分类与制表”（范围问题）
5. FRAGSTATS（我无法加载 shp 文件）
6. 光栅 --> 提取 --> 剪裁器不起作用（环自相交）

我听说 ArcMap 中的 Tabulate Area 可以做到这一点，但如果有一个开源解决方案，那就太好了。

Answer 1

我已经成功地用Python“rasterio”和“geopandas”做到了

它现在创建一个表，如下所示：

因为我没有找到类似 R“raster”中的 extract 命令的东西，所以它只需要 2 行，但不是计算半个晚上，现在一年只需要 2 分钟。 结果是一样的。它基于“https://gis.stackexchange.com/questions/260304/extract-raster-values-within-shapefile-with-pygeoprocessing-or-gdal/260380”中的“基因”思想

import rasterio
from rasterio.mask import mask
import geopandas as gpd
import pandas as pd

print('1. Read shapefile')

shape_fn = "D:/path/path/multypoly.shp"
raster_fn = "D:/path/path/class_1992.tif"

# set max and min class
raster_min = 10
raster_max = 230

output_dir = 'C:/Temp/'

write_zero_frequencies = True
show_plot = False

shapefile = gpd.read_file(shape_fn)

# extract the geometries in GeoJSON format
geoms = shapefile.geometry.values # list of shapely geometries
records = shapefile.values

with rasterio.open(raster_fn) as src:

     print('nodata value:', src.nodata)

     idx_area = 0
     # for upslope_area in geoms:
     for index, row in shapefile.iterrows():

          upslope_area = row['geometry']
          lake_id = row['ABC_ID']
          print('\n', idx_area, lake_id, '\n')

          # transform to GeJSON format
          from shapely.geometry import mapping
          mapped_geom = [mapping(upslope_area)]

          print('2. Cropping raster values')
          # extract the raster values values within the polygon
          out_image, out_transform = mask(src, mapped_geom, crop=True)

          # no data values of the original raster
          no_data=src.nodata

          # extract the values of the masked array
          data = out_image.data[0]

          # extract the row, columns of the valid values
          import numpy as np
          # row, col = np.where(data != no_data)
          clas = np.extract(data != no_data, data)

          # from rasterio import Affine # or from affine import Affine
          # T1 = out_transform * Affine.translation(0.5, 0.5) # reference the pixel centre
          # rc2xy = lambda r, c: (c, r) * T1

          # d = gpd.GeoDataFrame({'col':col,'row':row,'clas':clas})

          range_min = raster_min    # min(clas)
          range_max = raster_max    # max(clas)

          classes = range(range_min, range_max + 2)

          frequencies, class_limits = np.histogram(clas,
                                                   bins=classes,
                                                   range=[range_min, range_max])

          if idx_area == 0:
               # data_frame = gpd.GeoDataFrame({'freq_' + str(lake_id):frequencies})
               data_frame = pd.DataFrame({'freq_' + str(lake_id): frequencies})
               data_frame.index = class_limits[:-1]
          else:
               data_frame['freq_' + str(lake_id)] = frequencies
      idx_area += 1

 print(data_frame)
 data_frame.to_csv(output_dir + 'upslope_area_1992.csv', sep='\t')

Answer 2

我可以使用 R 命令提取来完成此操作，并用表格对其进行总结，如“Spacedman”所解释的，请参阅：https://gis.stackexchange.com/questions/23614/get-raster-values-from-a-polygon-开源地理信息解决方案中的叠加

shapes <- readOGR("C://data/.../shape)
LClass_1992 <- raster("C://.../LClass_1992.tif")
value_list <- extract (LClass, shapes )

stats  <- lapply(value_list,table)
[[354]]

10   11   30   40   60   70   80   90  100  110  130  150  180  190  200  201  210
67  303  233  450 1021 8241   65 6461 2823   88 6396    5   35  125   80   70 1027

但是需要很长时间（半夜）。 我会尝试用Python来做，也许会更快。 也许有人做过类似的事情并且可以分享代码。