计算多个重叠多边形（形状文件）内的栅格景观比例（百分比）？

Question

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

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

光栅包含类从 10 - 220 多边形有 ABC_ID 从 1 - 449

一年看起来像： classification and shape example

我需要创建一个 table 比如：

example table

我已经尝试通过以下方法实现此目的：

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

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

Answer 1

我可以使用 R 命令提取并使用 table 进行总结，如 "Spacedman" 所解释的，请参阅：https://gis.stackexchange.com/questions/23614/get-raster-values-from-a-polygon-overlay-in-opensource-gis-solutions

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 来做，也许会更快。也许有人做过类似的事情，可以分享代码。

Answer 2

我已经用 Python "rasterio" 和 "geopandas"

做到了

它现在创建一个 table 像： example result

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

的 "gene" 的想法

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')

计算多个重叠多边形（形状文件）内的栅格景观比例（百分比）？

Calculate raster landscape proportion (percentage) within multiple overlaping polygon (shapefiles)?

statistics

classification

polygon

raster

qgis