如何利用和操作“RasterBrick”中的土地覆盖图以进行 R 中的面积加权统计？答案

【问题标题】：How to utilize and manipulate land coverage map in `RasterBrick` for area-weighted statistics in R?如何利用和操作“RasterBrick”中的土地覆盖图以进行 R 中的面积加权统计？
【发布时间】：2018-05-09 09:13:55
【问题描述】：

我有TIF 格式的土地覆盖图，大概用于计算德国的区域加权年平均温度。我从这里 (direct download link of land coverage map for Europe) 下载了这个土地覆盖地图数据。特别是，我打算提取城市、农业区的土地/土壤覆盖数据，反之亦然。在我的第一步中，我使用raster 包导入了这些土地覆盖数据。下面是我的 R 脚本：

library(raster)
library(R.utils)
library(maptools)

url = "https://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/LUISA/PrimaryOutput/Europe/REF-2014/JRC_LUISA_Input_Corine_land_cover_2006_r_ref_2014.zip"
download.file(url, basename(url))
gunzip(basename(url))
rname <- list.files(getwd(), "tif$")

land_cover = raster::brick("~/LUISA_CLC_land_coverage/clc06_r.tif")

到目前为止，我可以在 R 中的RasterBrick 对象中导入原始土地覆盖图。请注意，原始地图覆盖了整个欧洲，所以我必须裁剪我只感兴趣的区域。为此，我使用maptools 和raster 包来裁剪地图。下面是 R 脚本：

data(wrld_simpl)
germany <- wrld_simpl[wrld_simpl@data$NAME == "Germany",]
germany <- spTransform(germany, CRSobj = land_cover@crs)
germany_land <- crop(land_cover, germany)

但是，我假设RasterBrick 对象中的这张裁剪后的土地覆盖图最好位于网格shapefile 中，具有非常高的分辨率，但它是如何实现的？任何的想法？

提出这个问题的要点是我需要检索城市、农业区的所有土地/土壤覆盖数据，并将这些信息与相应的德国 NUTS-s 级别 shapefile (download link of Germany level 3 shapefile) 匹配。

我真的不知道如何利用这张土地覆盖地图中的数据来计算面积加权的年平均温度。或许，一种可能的方法是检索城市的土地/土壤覆盖率、农业区数据，然后从德国 NUTS-3 级别的 shapefile 中找到匹配项。

这里是如何获取德国的 NUTS-3 shapefile（R 脚本如何在 R 中获取德国的 NUTS-3 地区的 shapefile）：

library(maptools)
library(rgdal)
library(R.utils)

url = "http://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2013-03m.shp.zip"
download.file(url, basename(url))
gunzip(basename(url))

getwd()
setwd("~/ref-nuts-2013-03m.shp/")
list.files(pattern = 'NUTS_RG_03M_2013.*.shp$')

eu <- readOGR(dsn = getwd(), layer ="NUTS_RG_03M_2013_4326_LEVL_2")
Germany_NUTS3 <- eu[eu@data$CNTR_CODE == "DE",]

所以使用这个 Gernamny 的 NUTS-3 shapefile Germany_NUTS3，我想提取城市、农业区的所有土地/土壤覆盖数据，反之亦然。

如果从RasterBrick 的土地覆盖图中提取数据，我如何在 R 中实现这一点？那可行吗？任何有效的解决方法来完成这项工作？任何的想法？

【问题讨论】：

什么样的信息？土地覆盖数据只是一个为每个像素分配一个类的栅格 - 所以从我的脑海中，您可以提取 shapefile 中每个要素的像素并计算每个类的百分比
最终，这将取决于您的目标是什么……我只是假设您有 NUTS 区域和土地覆盖栅格，这将是一种将信息汇总到每个区域的自然方式。因此，假设其中一个区域被 1000 个像素覆盖，如果其中 100 个被标记为“裸土”，您可以争辩说该区域的 10% 被裸土覆盖（当然，这个结论将取决于您的质量土地覆盖数据，像素越多，结论越稳健）
在光栅中，每个像素都有一个不同的类......因此您可以提取每个多边形的所有像素（查看raster 包的extract 函数）。这可以让您了解每个多边形内土地覆盖类别的分布。此外，您可以通过像素计数等方法来估计面积，这与光栅的分辨率直接相关……但这正在变成理论讨论，而 SO 并不是真正适合这样做的地方。我建议你阅读一些文献和谷歌搜索（例如从土地覆盖分类估计土地覆盖面积）
Like this or this ... 如果你对自己想要做什么有一个清晰的想法，你可以考虑程序化的实现。希望有帮助！
让我们continue this discussion in chat。

标签： r geospatial raster data-manipulation

【解决方案1】：

正如我们在 cmets 和聊天中所讨论的，这将是一种快速而肮脏的方法（其中 JRC 方法肯定是一种更好的方法，但也需要更多的努力）

所以你有你的土地覆盖land_cover 和你在德国的 NUTS 区域Germany_NUTS3：

# you can take the raster function since it's only one layer

land_cover <- raster::raster("~/LUISA_CLC_land_coverage/clc06_r.tif")

eu <- readOGR(dsn = getwd(), layer ="NUTS_RG_03M_2013_4326_LEVL_2")
Germany_NUTS3 <- eu[eu@data$CNTR_CODE == "DE",]

# you can use Germany_NUTS straight for cropping the landcover, so we'll project it to the raster's coordinate system

Germany_NUTS3_projected <- spTransform(Germany_NUTS3, CRSobj = land_cover@crs)

land_cover_germany <- crop(land_cover, Germany_NUTS3_projected)

现在您可以使用光栅包中的extract 提取 NUTS 区域的土地覆盖像素。

注意：这可能需要一些时间，尤其是在区域或栅格很大或者您有很多多边形的情况下。如果您需要重复执行此操作，您可能需要使用不同的方法。

例如，我将针对 NUTS 地区之一进行此操作：

pixel_extract <- raster::extract(land_cover_germany,Germany_NUTS3_projected[2,])

如果您一次使用多个多边形，pixel_extract 将是一个带有像素值的向量列表，每个元素代表一个不同的多边形。

在本示例中，只有一个元素显示了该多边形内像素的土地覆盖类值：

> head(pixel_extract)
[[1]]
   [1] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
  [24] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
  [47] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
  [70] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
  [93] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
 [116] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 23 23 24 24 24 24 24 24 24
 [139] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24  4 ...

现在，要导出您感兴趣的类所覆盖的区域，您需要计算像素，然后将它们与单个像素的面积相乘。由于单个像素的分辨率为 100 x 100 米，因此面积为 10000 平方米。

要确定所需类别的土地覆盖值，我们可以查看栅格随附的LUISA_legend.xls 文件：

   GRID_CODE CLC_CODE LABEL
    1   111 Artificial surfaces
    2   112 Artificial surfaces
    3   121 Artificial surfaces
    4   122 Artificial surfaces
    5   123 Artificial surfaces
    6   124 Artificial surfaces
    7   131 Artificial surfaces
    8   132 Artificial surfaces
    9   133 Artificial surfaces
    10  141 Artificial surfaces
    11  142 Artificial surfaces
    12  211 Agricultural areas
    13  212 Agricultural areas
    14  213 Agricultural areas
    15  221 Agricultural areas
    16  222 Agricultural areas
    17  223 Agricultural areas
    18  231 Agricultural areas
    19  241 Agricultural areas
    20  242 Agricultural areas
    21  243 Agricultural areas
    22  244 Agricultural areas

所以要计算像素，我们只需要查看人造表面的值在 1 到 11 之间，对于农业表面的值在 12 到 22 之间。要获得面积“估计”，我们可以用单个像素的面积计算像素数。

areas <-data.frame(ARTIFICIAL_AREA=sum(pixel_extract[[1]]>=1 &  pixel_extract[[1]]<=11) * (100*100),
           AGRICULTURE_AREA=sum(pixel_extract[[1]]>=12 &  pixel_extract[[1]]<=22) * (100*100))

这会为您提供以平方米为单位的面积估算值：

> areas
  ARTIFICIAL_AREA AGRICULTURE_AREA
1       954030000       9282970000

如果pixel_extract 是一个每个多边形都有一个元素的列表（也就是说，如果您使用了完整的 shapefile），您可以使用 lappy 计算面积并使用 do.call(rbind, 创建单个数据框：

areas <- lapply(pixel_extract, function(x) data.frame(ARTIFICIAL_AREA=sum(x >=1 &  x <=11) * (100*100),
                   AGRICULTURE_AREA=sum(x>=12 & x<=22) * (100*100))

areas_df <- do.call(rbind,areas)

【讨论】：