如何使用 raster 包中的函数 crop(x, y, ...) 以便输出是关于 x (shapefile, extent object) 中包含的 y (raster) 的信息

Question

我对 R 有点陌生，尤其是在 R 中使用 GIS，所以我希望我能正确解释这一点。

我想从包 raster 中获取函数 crop (x,y...) 以合并/覆盖（不确定要使用的正确表达式是什么）带有 shapefile 的光栅文件。本质上，shapefile 是瑞典的 5 x 5 公里网格，栅格是瑞典的土地利用栅格。通过使用裁剪函数，我想得到一个结果，对于 shapefile 中的每个 5x5 平方千米，提供从栅格文件中提取的有关每个正方形中土地利用类型的信息。

但是，当我使用crop，然后从我裁剪的“crop (landuse, ekrut)”中写入.csv（参见下面的代码）时，我基本上只是将 ekrut（shapefile）作为输出 csv 文件，没有列添加了土地利用栅格，指示哪个广场具有哪种土地利用类别。

很抱歉，我不能在这里分享实际的 shapefile，如果它会有所作为，可以从这里下载土地利用数据： http://gpt.vic-metria.nu/data/land/NMD/NMD2018_basskikt_ogeneraliserad_Sverige_v1_0.zip

这是我到目前为止尝试做的事情

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

这是 GIS 文件的坐标系/投影。每个坐标 系统有一个 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 文件，没有添加土地利用栅格的列，指示哪个广场有哪个土地利用类。

我想要一个 .csv，对于每个广场（其中有 19 191 个），我可以提供有关该广场存在何种土地利用类型的信息。

如果有更好的方法来解决这个问题，请指出。我希望我的解释很清楚！

谢谢。

Answer 1

即使您的代码正常工作，您仍然无法将此输出导出到 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 = '')

Answer 2

您可以提取每个多边形的所有类，为每个多边形制作一个表格，然后将它们组合起来。

以 Majid 的示例数据为基础。

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)