I'm trying to find out the land use type for a set of coordinates that define the location of plant species across Europe. However I'm stuck in the process of assigning the land use to the respective coordinates. Any advice would be more than welcome!
First, I download the land use raster file from here: 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")
Then, I read the Corine land use classes and rename the levels of the raster file with this classes
#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)
This is a small subset of coordinates from my full list of coordinates (over 200.000 in total):
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)
But now I really do not know how to proceed and create the final dataframe with land use type given the matching with the raster file
My intention is to add a 4th column as follows after the matching of my coordinates with the land use type of the raster file.
#Example of how this fourth column would be like:
d$land_use <- c("Olive groves", "Olive groves", "Vineyards", "Pastures")
The data (another file from the same website).
library(terra)
r <- rast("U2018_CLC2018_V2020_20u1.tif")
As you can see, r knows about the class labels.
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
Here are some example points for use with 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
Or with your lon/lat points (you had the names reversed!), first transform these to the coordinate reference system of the land use raster:
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
If you want the land use code instead
activeCat(r) <- "CODE_18"
extract(r, vv)
# ID CODE_18
#1 1 242
#2 2 122
#3 3 231
#4 4 121