Как использовать функцию crop(x, y, ...) из растрового пакета, чтобы на выходе была информация о y (растре), содержащаяся в x (шейп-файл, объект экстента)
Я немного новичок в R и особенно в работе с GIS в R, поэтому я надеюсь, что объясню это правильно.
Я хочу получить функцию обрезки (x,y...) из растрового пакета для слияния / наложения (не знаю, какое правильное выражение использовать) растровый файл с шейп-файлом. По сути, шейп-файл представляет собой сетку Швеции размером 5 × 5 км, а растр - это растровый растр Швеции. Используя функцию обрезки, я хочу получить результат, который для каждого квадрата 5 x 5 км из шейп-файла дает информацию, извлеченную из растрового файла, о типе землепользования в каждом квадрате.
Тем не менее, когда я использую crop, а затем записываю.csv из того, что я обрезал "crop (landuse, ekrut)" (см. Код ниже), я в основном просто получаю ekrut (shapefile) в качестве выходного файла csv, столбцы из растр землепользования, указывающий, какой квадрат имеет какой класс землепользования.
Извините, но я не могу поделиться с вами фактическим шейп-файлом, если это как-то изменится, данные о землепользовании можно скачать здесь: 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 из того, что я обрезал, я в основном просто получаю ekrut (shapefile) в качестве выходного csv-файла, не добавляются столбцы растра landuse, указывающие, какой квадрат имеет какой класс землепользования.
Я хотел бы иметь.csv, который для каждого квадрата (их 19 191) дает мне информацию о том, какой тип землепользования присутствует в этом квадрате.
Если есть лучший способ приблизиться к этому, пожалуйста, укажите. Я надеюсь, что я был ясен с моим объяснением!
Благодарю.
2 ответа
Вы можете извлечь все классы для каждого многоугольника, создать таблицу для каждого многоугольника, а затем объединить их.
Основываясь на данных примера Маджида.
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)
Даже если ваш код работал, вы все равно не можете экспортировать этот вывод в 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)
ht tps:https://stackru.com/images/c88bb69b84fcc6dcddf71ac65930219607c00bb7.jpg
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 = '')