Как использовать и манипулировать картой покрытия земли в RasterBrick для взвешенной по площади статистики в R?

Question

У меня была карта покрытия земли в формате TIF, которая предположительно использовалась для расчета средневзвешенной по площади среднегодовой температуры для Германии. Я скачал данные карты покрытия земли здесь ( прямая ссылка для скачивания карты покрытия земли для Европы ). В частности, я намерен извлечь данные о земельном / почвенном покрове для города, сельскохозяйственной территории и наоборот. На первом этапе я импортировал эти данные о земельном покрытии с пакетом 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")

пока что я могу импортировать исходную карту покрытия земли в объекте RasterBrick в R. Обратите внимание, что исходная карта покрыла всю Европу, поэтому я Я должен обрезать регионы, которые меня интересуют. Для этого я использовал maptools и raster package для отсечения карты. Вот сценарий 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 Германии ( ссылка для скачивания уровня Германии 3 шейп файла ).

У меня действительно нет четкого представления о том, как использовать данные на этой карте земельного покрытия для расчета среднегодовой температуры, взвешенной по площади. Возможно, возможный подход мог бы заключаться в том, чтобы получить земельный / почвенный охват города, данные по сельскохозяйственным районам, а затем найти соответствие из шейп-файла уровня NUTS-3 Германии.

Вот как получить шейп-файл NUTS-3 для Германии (скрипт R, как получить шейп-файл для регионов Германии NUTS-3 в R):

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, 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, используя extract из растрового пакета.

ВНИМАНИЕ: это может занять некоторое время, особенно если область или растр велики или у вас много полигонов. Если вам нужно делать это повторно, вы можете использовать другой подход.

В качестве примера я сделаю это для одного из регионов NUTS:

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

Если вы используете больше полигонов одновременно, pixel_extract будет списком векторов со значениями пикселей, причем каждый элемент представляет отдельный полигон.

В этом примерном случае есть только один элемент, показывающий значения класса Landcover для пикселей в этом многоугольнике:

> 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 м2.

Чтобы определить значения земного покрова для желаемых классов, мы можем взглянуть на файл 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 является списком с элементом на полигон (то есть, если вы использовали полный шейп-файл), вы можете вычислить области с помощью 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)