Как использовать функцию crop (x, y, ...) из растрового пакета, чтобы на выходе была информация о y (растре), содержащаяся в x (шейп-файл, объект экстента)

Question

Я немного новичок в R и особенно в работе с ГИС в R, поэтому я надеюсь, что объясню это правильно.

Я хочу получить функцию обрезки (x, y ...) из растрового пакета для слияния / наложения (не знаю, какое правильное выражение использовать) растровый файл с шейп-файлом.По сути, шейп-файл представляет собой сетку Швеции размером 5 × 5 км, а растр - это растровый растр Швеции.Используя функцию обрезки, я хочу получить результат, который для каждого квадрата 5x5 км из шейп-файла дает извлеченную из растрового файла информацию о типе землепользования в каждом квадрате.

Однако, когда я использую обрезкуи затем напишите .csv из того, что я обрезал 'crop (landuse, ekrut) "(см. код ниже), я в основном просто получаю ekrut (shapefile) в качестве выходного csv-файла, нет добавленных столбцов растра landuse, указывающего, какой квадрат имееткакой класс землепользования.

Извините, но я не могу поделиться фактическим шейп-файлом здесь, если это как-то изменится, данные о землепользовании можно скачать здесь: http://gpt.vic -metria.nu/data/land/NMD/NMD2018_basskikt_ogeneraliserad_Sverige_v1_0.zip

Вот что я пытался сделать до сих пор

library (raster)
library (rgdal)
library (sp)

это система координат / проекция для файла ГИС. Каждая координата.система имеет код epsg (http://spatialreference.org/).

sweref.def <- "+init=epsg:3006" 

# read in shapefile
ekrut <- readOGR (dsn   = "//storage-al.slu.se/student$/nilc0001/Desktop/Nina/Ekrut", 
                  layer = "ekrut_5x5_flat", 
                  p4s   = sweref.def) 

ekrut
# class       : SpatialPolygonsDataFrame 
# features    : 19192 
# extent      : 266333, 921700, 6131565, 7675329  (xmin, xmax, ymin, ymax)
# coord. ref. : +init=epsg:3006 +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# variables   : 7
# names       :    AREA, PERIMETER,  GGD_, GGD_ID,     BK, Bk_num, BK_flat 
# min values  : 2.5e+07,     20000,     2,      0, 10A 0f, 012 77,   10A0f 
# max values  : 2.5e+07,     20000, 19208,  19230,  9J 9d, 329 48,    9J9d 

#read in raster
landuse <- raster ("nmd2018bas_ogeneraliserad_v1_0.tif")

landuse

# class       : RasterLayer 
# dimensions  : 157992, 71273, 11260563816  (nrow, ncol, ncell)
# resolution  : 10, 10  (x, y)
# extent      : 208450, 921180, 6091140, 7671060  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# data source : //storage- al.slu.se/student$/nilc0001/Desktop/Nina/landuse/nmd2018bas_ogeneraliserad_v1_0.tif 
# names       : nmd2018bas_ogeneraliserad_v1_0 
# values      : 0, 128  (min, max)
# attributes  :
#         ID      COUNT Opacity Klass
#  from:   0 5204803484       0      
#  to  : 255          0       0   

#first few rows of attribute table of landuse
levels (landuse)   

# [[1]]
#      ID      COUNT Opacity                                              Klass
# 1     0 5204803484       0                                                   
# 2     1          0     255                                                   
# 3     2  382320369     255                                      Öppen våtmark
# 4     3  258249590     255                                           Åkermark

# crop and write csv
landuse.ekrut <- crop (landuse, ekrut)
write.csv (landuse.ekrut,"landuse.ek.csv")

Как я уже говорил, когда я использую кадрирование и затем пишу .csv из того, что я обрезал, я в основном просто получаюt ekrut (shapefile) в качестве выходного CSV-файла, нет добавленных столбцов растра землепользования, указывающих, какой квадрат имеет какой класс землепользования.

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

Если есть лучший способ приблизиться к этому, пожалуйста, укажите.Надеюсь, я объяснил свое объяснение!

Спасибо.

Answer 1

Вы можете извлечь все классы для каждого многоугольника, составить таблицу для каждого многоугольника, а затем объединить их.

Опираясь на пример данных Маджида.

library(raster)
set.seed(1523)
r <- raster(ncol=9, nrow=6, vals=sample(10, 9*6, replace=TRUE))
rr <- aggregate(raster(r), fact=3)
pols <- rasterToPolygons(rr) # transform into a polygon
plot(r); lines(pols);  text(pols, 1:length(pols))

Извлечение, создание таблицы и объединение

e <- extract(r, pols)
x <- lapply(e, table)
y <- lapply(1:length(x), function(i) data.frame(pol=i, as.data.frame(x[[i]])))
z <- do.call(rbind, y)

head(z, 10)
#   pol Var1 Freq
#1    1    1    2
#2    1    2    2
#3    1    3    2
#4    1    6    2
#5    1    9    1
#6    2    1    1
#7    2    2    1
#8    2    4    2
#9    2    6    2
#10   2    7    1

и сохранить в CSV

write.csv(z, "test.csv", row.names=FALSE)

Answer 2

Даже если ваш код работал, вы все равно не сможете экспортировать этот вывод в CSV, если не используете что-то вроде majority фильтра. Это в основном из-за того, что один полигон может содержать более одного землепользования.

Один из способов сделать это, как показано ниже:

library(raster)
library(stats)

#set.seed to have a reproducible example
set.seed(1523)

#making a raster layer
r <- raster(ncol=36, nrow=18, vals=floor(runif(648, 5.0, 7.5)))
rr <- aggregate(r, fact=3)
f.net <- rasterToPolygons(rr) # transform into a polygon

#simple ovelaying plot for visualisation
plot(r, main="Raster and the overlaying poly")
plot(f.net, col=rgb(1, 0, 0, alpha = 0.3), add=TRUE)

ext.r <- extract(r, f.net) # returns values per ploy

# a simple mode filter
mode <- function(x){ 
  ta = table(x)
  tam = max(ta)
  mod = as.numeric(names(ta)[ta == tam])
  mod = mod[1] 
  return(mod)
}

# each poly contains several values
head(ext.r)
# [[1]]
# [1] 6 5 5 6 5 6 5 6 5
# 
# [[2]]
# [1] 5 5 6 7 6 7 6 5 7
# 
# [[3]]
# [1] 5 5 7 6 6 7 6 6 7
# 
# [[4]]
# [1] 6 7 5 5 7 6 5 6 5
# 
# [[5]]
# [1] 5 6 5 5 5 5 6 5 5
# 
# [[6]]
# [1] 7 5 6 6 6 5 5 6 5


# a mode filter can be applied to pas majority values to each poly
mode.ext.r <- lapply(ext.r, mode)
head(mode.ext.r)
# [[1]]
# [1] 5
# 
# [[2]]
# [1] 5
# 
# [[3]]
# [1] 6
# 
# [[4]]
# [1] 5
# 
# [[5]]
# [1] 5
# 
# [[6]]
# [1] 5
write.table(mode.ext.r, file = 'mode_ext_r.csv', row.names = FALSE, na = '')