0

I would like to extract values from Corine Land Cover (CLC) 2018 raster.

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

r <- raster('U2018_CLC2018_V2020_20u1.tif')
#Warning message:
#In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#  Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition

I got a Warning message, probably due to a lack of a +datum parameter in the proj4 string. The Coordinate Reference System (CRS) for the CLC 2018 is EPSG:3035 (ETRS89, LAEA). The crs from the raster 'r' is :

crs(r)
#CRS arguments:
# +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 

Then I load my 'data' file that is a 'data.frame' with 419389 rows and 15 columns including the 'Latitude' and 'Longitude'. I then run these few lines

lats <- c(data$Latitude)
lons <- c(data$Longitude) 
coords <- data.frame(x=lons,y=lats)
points <- SpatialPoints(coords, proj4string = r@crs) 
points <- na.omit(object = points)
values <- extract(r, points)

points
#class       : SpatialPoints 
#features    : 419389 
#extent      : -8.06205, 2.39726, 49.18919, 60.59815  (xmin, xmax, ymin, ymax)
#crs         : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

However, the 'values' return NAs.

I guess the problem is the projection, so I tried to reproject the raster but it doesn't work.

projection(r) <- crs("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +datum=WGS84 +units=m +no_defs")
1

You are mixing up the difference between setting a crs (declaration) of a spatial dataset and creating a new spatial dataset with a different crs (transformation, reprojection).

You need to correctly declare the crs of your points. You set them to the crs of r, but that is wrong --- at least the variable names suggest it should instead be longitude/latitude, so you can do:

points <- SpatialPoints(coords, proj4string =CRS("+proj=longlat +datum=WGS84"))

Since this crs does not match that of the raster data, you need to transform the points (transforming point data is fast and does not degenerate the data; that is not true for transforming raster data)

pp <- spTransform(points, "EPSG:3035")
 

And then you should be able to use extract.

1
  • Thank you very much! It works, and I got where I was wrong. Apr 5 '21 at 10:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.