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

Question

我有一张TIF格式的土地覆盖图，大概用来计算德国地区加权年平均温度。我从这里下载了这张土地覆盖图数据（欧洲土地覆盖图的直接下载链接）。特别是，我打算提取城市和农业区的土地/土壤覆盖率数据，反之亦然。在我的第一步中，我使用光栅软件包导入了这个土地覆盖率数据。下面是我的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级形状文件（德国3级形状文件的下载链接）相匹配。

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

以下是如何获取德国“NUTS-3”形状文件（R脚本如何在R中获取德国“NUTS-3 regions”形状文件）：

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中实现这一点？这可行吗？有什么有效的解决方法来完成这个任务吗？有什么想法吗？

Answer 1

正如我们在评论和聊天中所讨论的，这将是一种快速而肮脏的方法（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)

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

注意：这可能需要一些时间，尤其是如果区域或光栅很大或者你有很多多边形。如果你需要重复这样做，你可能想使用不同的方法。

例如，我将对其中一个坚果区域执行此操作：

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

如果一次使用多个多边形，“pixel\u 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×100米，因此面积为10000平方米。

要确定所需类别的土地覆盖值，我们可以查看LUISA\u图例。光栅随附的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)