Любой эффективный способ получить оценку площади с подходом подсчета пикселей для растра земного покрова в R?

Question

Я хочу получить оценку площади для каждого полигона, используя метод подсчета пикселей.Первоначальная карта земного покрова представляла собой растровые данные, которые снабжены отличительным классом, назначенным для каждого пикселя ( легенда карты земного покрова ).Однако, чтобы получить оценку площади с помощью подсчета пикселей, мне не интуитивно понятно.Возможно, первое, что я могу сделать, это извлечь все пиксели для каждого многоугольника, которые представляют информацию о распределении классов земного покрова в каждом многоугольнике.После этого мне нужно использовать метод подсчета пикселей, чтобы получить оценку площади и агрегировать весь земельный / почвенный охват города, сельскохозяйственную площадь для каждого полигона.Для меня использование подсчета пикселей для получения оценки площади не является простым и не имеет четкого представления о том, как это сделать в R.

Обратите внимание, что оригинальная карта земельного покрова довольно большая (около 92 mb)и трудно привести воспроизводимый пример для растрового покрова, так что прости меня за такие неудобства.Вот скрипт R для получения растра почвенно-растительного покрова:

library(raster)
library(R.utils)

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::raster("~/LUISA_CLC_land_coverage/clc06_r.tif")

Я хочу объединить всю информацию о земельном покрытии / почве для каждого полигона (всего 403 полигона, которые будут агрегировать информацию о земном покрове);Вот полигоны, которые я собираюсь использовать для оценки площади: шейп-файл с полигонами доступен на лету :

Я обрезал оригиналрастр Landcover следующим образом:

deu_shp = maptools::readShapePoly("~/adm2.shp",
                                  proj4string=crswgs84,verbose=TRUE)
deu_proj <- spTransform(deu_shp, CRSobj = land_cover@crs)
land_cover_deu <- crop(land_cover, deu_proj )

Чтобы понять оценку площади с подсчетом пикселей, я погуглил связанную статью и нашел это: Использование дистанционного зондирования для оценки площади посева и идеи, представленной в этой статьеблизко соответствует моим интересам, но реализовать представленный метод является чисто теоретическим, и мне сложно это сделать в R.Я вполне в порядке, если есть какой-либо быстрый и грязный способ получить оценку площади с подсчетом пикселей, где информация о земельном покрытии / почве (например, город, лес, сельскохозяйственная земля) должна быть агрегирована для каждого многоугольника (доступен шейп-файл с многоугольниками)на лету ).

для выделения пикселей, я мог бы просто использовать raster::extract следующим образом (приведенный ниже код является пробным):

lapply(deu_proj, function(x) {
    pixel_extract = raster::extract(land_cover_deu, deu_proj[x,])
    pixel_extract= as.data.frame(pixel_extract)
})

и выше, простое извлечение пикселей не очень эффективно для исходного земного покроварастр.Я не знаю, как сделать подсчет пикселей для извлеченного пикселя в каждом многоугольнике и получить соответствующую оценку площади (например, площадь земли в городе, лес, сельское хозяйство и т. Д.).

Любую идею сделатьэто случилось в R?Как получить оценку площади с использованием метода подсчета пикселей для данного растра земного покрова?Возможно ли объединить информацию о земельном покрытии / почве для каждого полигона?Заранее спасибо

Answer 1

Вот рабочий процесс для получения количества пикселей каждого класса растительного покрова по полигонам (операция, называемая областью табуляции). Идея состоит в том, чтобы растеризовать шейп-файл, используя разрешение и размер растрового покрова. Затем вычислите область табулирования, используя очень эффективную функцию на основе data.table, которая подсчитывает количество пикселей каждого класса в каждом многоугольнике.

# add ID field to the shapefile
deu_proj@data$ID <- 1:nrow(deu_proj@data)

# extract extent and resolution of land cover raster
ext <- extent(land_cover_deu)
ext <- paste(ext[1], ext[3], ext[2], ext[4])
res <- paste(res(land_cover_deu)[1], res(land_cover_deu)[2])

# rasterize shapefile using gdal (more efficent than rasterize from raster package)
# you can change GDAL_CACHEMAX according to your RAM
system(paste0("gdal_rasterize --config GDAL_CACHEMAX 8000 -a ID -te ",
               ext," -tr ", res,
               " -ot Int16 /home/deu_proj.shp /home/zones.tif"))

# load raster
zones <- raster("/home/zones.tif")

# zonal stats using myZonal function
Zstat <- myZonal(land_cover_deu, zones)

# reshape output    
results <- data.table::dcast(Zstat, z ~ vals)
colnames(results)[1] <- "ID"

# merge data
deu_proj@data <- plyr::join(deu_proj@data, results, by="ID")

# show results
head(deu_proj@data[, c(8, 17:ncol(deu_proj@data))])

           NAME_2  0   1    2    3    4  5   6   7  8   9 10  11    12  15   16    18    20   21 22    23    24    25  26
1 Alb-Donau-Kreis NA 241 4504  781 1317 NA  NA 357 NA  NA NA 133 57460  NA   NA  7212 19899 1627 NA 16403  6628 15248 135
2       Böblingen NA 369 4947 1599 1140 NA  NA 149 87 116 98 438 17845  NA 1712  1717  4742 3079 NA  3988  2286 14557  NA
3     Baden-Baden NA  45  791  150  306 NA  83  33 NA  NA 38  41   695 338  602   721   468   92 NA  1090  2272  4401  NA
4        Biberach NA 201 4681  462 1264 NA 152 312 NA  48 NA 131 46163  NA   29 23780 19126  920 NA  2291 28679  5140  NA
5   Bodenseekreis NA 268 3219  561  814  7 169  81 NA  NA NA  76  9482 418 5403  4750 19105  804 NA    96  9151  8797  NA
6        Bodensee NA  NA   13    1   NA  9  NA  NA NA  NA NA   3    NA  NA   NA     1     2   NA NA    NA    NA     4  NA
  27   29 30 31 32 34   35  36 37 39 40  41 42 43   45 46   128
1 NA  581 NA NA NA NA  104  NA NA NA NA 397 NA NA 2643 40    NA
2 NA  801 NA NA NA NA   NA  NA NA NA NA  NA NA NA 2001 58    NA
3 NA 1117 NA NA NA NA  157  NA NA NA NA  79 NA NA  486  1    NA
4 NA 2063 NA NA NA NA 1105 273 NA NA NA 384 NA NA 3760 25    NA
5 NA  240 NA NA NA NA  299  NA NA NA NA 193 NA NA 2476 66    45
6 NA   NA NA NA NA NA    4  NA NA NA NA 415 NA NA   20  1 27563

Функция Zonal Stats (адаптирована с здесь )

myZonal <- function (x, z, digits = 0, na.rm = TRUE) { 
  vals <- raster::getValues(x) 
  zones <- round(raster::getValues(z), digits = digits) 
  rDT <- data.table::data.table(vals, z = zones) 
  plyr::count(rDT) }

Пример данных

# I have loaded the shapefile using getData but you can use your own shp.
# The only difference will be in the column numbers of "show results" step.
deu_shp <- getData("GADM", country="DEU", level=2)