如何将 Corine 的土地覆盖类型分配和匹配到具有一组经纬度坐标的数据框？

Question

我正在尝试找出一组坐标的土地利用类型，这些坐标定义了整个欧洲植物物种的位置。但是，我陷入了将土地使用分配给相应坐标的过程中。任何建议都将非常受欢迎！

首先，我从这里下载土地利用栅格文件：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)

但现在我真的不知道如何继续并创建与栅格文件匹配的土地利用类型的最终数据框

我的意图是在我的坐标与光栅文件的土地利用类型匹配后添加如下第 4 列。

#Example of how this fourth column would be like:
d$land_use <- c("Olive groves", "Olive groves", "Vineyards", "Pastures")

Answer 1

数据（同一网站的另一个文件）。

library(terra)
r <- rast("U2018_CLC2018_V2020_20u1.tif")

如您所见，r 知道 class 个标签。

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

或者用你的 lon/lat 点（你的名字颠倒了！），首先将它们转换为土地利用栅格的坐标参考系统：

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

如何将 Corine 的土地覆盖类型分配和匹配到具有一组经纬度坐标的数据框？

How to assign and match land cover type from Corine to a dataframe with a set of lon lat coordinates?

gis

r