我正在尝试找出一组定义欧洲植物物种位置的坐标的土地利用类型。然而,我陷入了将土地用途分配给相应坐标的过程中。任何建议都非常受欢迎!
首先,我从这里下载土地利用栅格文件:https://land.copernicus.eu/pan-european/corine-land-cover
#Read raster file (year 2006 but could be any)
clc <- raster("U2006_CLC2000_V2020_20u1.tif")
然后,我读取了 Corine 土地利用类并使用该类重命名了栅格文件的级别
#Read Corine classes
clc_classes <- foreign::read.dbf("CLC_1990/DATA/U2006_CLC2000_V2020_20u1.tif.vat.dbf",
as.is = TRUE) %>%dplyr::select(value = Value,landcov = LABEL3)
这是我的完整坐标列表(总共超过 200.000 个)中的一小部分坐标:
lon <- c("51.105", "51.195", "51.188", "51.239")
lat <- c("4.392", "4.395", "4.896", "4.468")
sp <- c("sp1","sp2", "sp3","sp4")
#Create minimal dataframe
d <- data.frame(lon,lat,sp)
但现在我真的不知道如何继续并根据与栅格文件的匹配创建具有土地利用类型的最终数据框
我的目的是在我的坐标与栅格文件的土地利用类型匹配后添加第四列,如下所示。
#Example of how this fourth column would be like:
d$land_use <- c("Olive groves", "Olive groves", "Vineyards", "Pastures")
数据(同一网站的另一个文件)。
library(terra)
r <- rast("U2018_CLC2018_V2020_20u1.tif")
如您所见,
r
知道类标签。
r
#class : SpatRaster
#dimensions : 46000, 65000, 1 (nrow, ncol, nlyr)
#resolution : 100, 100 (x, y)
#extent : 9e+05, 7400000, 9e+05, 5500000 (xmin, xmax, ymin, ymax)
#coord. ref. : ETRS_1989_LAEA (EPSG:3035)
#source : U2018_CLC2018_V2020_20u1.tif
#color table : 1
#categories : LABEL3, Red, Green, Blue, CODE_18
#name : LABEL3
#min value : Continuous urban fabric
#max value : NODATA
head(cats(r)[[1]])
# Value LABEL3 Red Green Blue CODE_18
#1 1 Continuous urban fabric 0.9019608 0.0000000 0.3019608 111
#2 2 Discontinuous urban fabric 1.0000000 0.0000000 0.0000000 112
#3 3 Industrial or commercial units 0.8000000 0.3019608 0.9490196 121
#4 4 Road and rail networks and associated land 0.8000000 0.0000000 0.0000000 122
#5 5 Port areas 0.9019608 0.8000000 0.8000000 123
#6 6 Airports 0.9019608 0.8000000 0.9019608 124
以下是一些与
extract
一起使用的示例点
pts <- matrix(c(3819069, 3777007, 3775822, 2267450, 2302403, 2331431), ncol=2)
extract(r, pts)
# LABEL3
#1 Sea and ocean
#2 Complex cultivation patterns
#3 Natural grasslands
或者使用您的经/纬度点(您的名称颠倒了!),首先将它们转换为土地利用栅格的坐标参考系:
lat <- c(51.105, 51.195, 51.188, 51.239)
lon <- c(4.392, 4.395, 4.896, 4.468)
xy <- cbind(lon, lat)
v <- vect(xy, crs="+proj=longlat")
vv <- project(v, crs(r))
extract(r, vv)
# ID LABEL3
#1 1 Complex cultivation patterns
#2 2 Road and rail networks and associated land
#3 3 Pastures
#4 4 Industrial or commercial units
如果您想要土地使用代码
activeCat(r) <- "CODE_18"
extract(r, vv)
# ID CODE_18
#1 1 242
#2 2 122
#3 3 231
#4 4 121
罗伯特的回答对于像我这样没有线索的人来说真的很棒。 然而,当我尝试实现同样的目标时,我得到了意想不到的结果。
当我跑步时
library(terra)
data <- rast("2018_CLC2018_V2020_20u1.tif")
data
我没有标签:
class : SpatRaster
dimensions : 46000, 65000, 1 (nrow, ncol, nlyr)
resolution : 100, 100 (x, y)
extent : 9e+05, 7400000, 9e+05, 5500000 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS_1989_LAEA (EPSG:3035)
source : U2018_CLC2018_V2020_20u1.tif
name : U2018_CLC2018_V2020_20u1
min value : 1
max value : 48
因此提取的结果在 1 - 48 范围内。 这是我在网站上找到的唯一与命名匹配的文件。 我很抱歉将此添加为答案,但我无法发表评论。
嗨,我正在尝试复制此内容,尽管我的 .tif 似乎没有正确加载,因为我收到此错误消息:
class : SpatRaster
dimensions : 5972, 5020, 1 (nrow, ncol, nlyr)
resolution : 100, 100 (x, y)
extent : 2863300, 3365300, 3211800, 3809000 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
source : U2018_CLC2018_V2020_20u1.tif
name : U2018_CLC2018_V2020_20u1
Warning message:
In class(object) <- "environment" :
Setting class(x) to "environment" sets attribute to NULL; result will no longer be an S4 object
如有任何建议,我们将不胜感激!