Extracting values with SpatialPointsDataFrame from Raster returns only NAs

Question

I am a beginner in using R for spatial data analysis and want to simply extract some values for certain locations (with lon/lat positions) from a raster stack. However, when using extract(), it returns only NAs, so I guess I made a mistake with the projection.

I loaded the lon/lat-locations into a SpatialPointsDataFrame, changed the lon/lat into numeric and set the crs to EPSG:4326. The raster stack comes with the crs EPSG:3035, so I used "sftransform()" to project the lon/lat-locations accordingly. "compareCRS" for the two datasets gave the output "TRUE", so the projection supposedly worked. When I try to use extract(), however, it returns only NA's.

When I put the coordinates into Google Maps, they are mapped correctly, with mapview() in R, however, they seem to be completely off grid, so I supposed something about the projection went wrong, but I have no idea what. I hope that somebody here can help me out!

Loading in the raster stack named "texture"

> print(texture)
class      : RasterStack 
dimensions : 8073, 7781, 62816013, 3  (nrow, ncol, ncell, nlayers)
resolution : 500, 500  (x, y)
extent     : 2635700, 6526200, 1385700, 5422200  (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 
names      : Clay_eu23, Sand_eu23, Silt_eu23 
min values :         0,         0,         0 
max values :  76.47613, 100.00000,  91.61756

Loading in the DF with Lon/Lat Locations

  Datalogger        Lat       Lon
1        DL3 48.5932246 9.5576115
2        DL4 49.0160648 9.1887693
3        DL2 48.7100801 9.2141318
4       DL10 49.0038758 8.4934965
5        DL5 49.0034701 8.4954198
6        DL6 48.4183515 8.8958118

Changing the Lon/Lat columns to numeric

locations$Lon <- as.numeric(locations$Lon)
locations$Lat <- as.numeric(locations$Lat)

Setting the crs for the coordinates and creating SpatialPointsDF

coordinates(locations) <- ~Lon+Lat
proj4string(locations) <- CRS("+init=epsg:4326")

> print(locations)
class       : SpatialPointsDataFrame 
features    : 10 
extent      : 1, 10, 1, 10  (xmin, xmax, ymin, ymax)
crs         : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : Datalogger 
min values  :        DL1 
max values  :        DL9

Reprojecting the locations to the raster stack's (texture's) crs

loc_proj <- spTransform(locations, proj4string(texture))

Checking if crs is correct

compareCRS(texture, loc_proj)

Trying to extract values from rasterstack

values <- raster::extract(texture, loc_proj)

For values I only get NAs for each of the three rasters in the stack, even though when I call summary(texture) it returns valid data, and also the locations as such are valid. What did I do wrong with the projections? Thank you very much already!!

Robert Hijmans Robert Hijmans · Accepted Answer · 2019-08-21T16:18:55

You did

locations$Lon <- as.numeric(locations$Lon)
locations$Lat <- as.numeric(locations$Lat)

You got

print(locations)
#class       : SpatialPointsDataFrame 
#features    : 10 
#extent      : 1, 10, 1, 10  (xmin, xmax, ymin, ymax)

That is, your longitude and latitude are between 1 and 10. This suggest that lon and lat were factor variables (for whatever reason); and that you made a mistake in changing them to numeric.

Consider this

x <- as.factor(c(9.5576115, 9.1887693, 9.2141318, 8.4934965, 8.4954198, 8.8958118))
as.numeric(x)
#[1] 6 4 5 1 2 3

as.numeric does not use the factor labels in the conversion. It cannot do that because these labels may contain letters. Instead, it gives you the indices

So, instead, you should do

as.numeric(as.character(x))
#[1] 9.557612 9.188769 9.214132 8.493496 8.495420 8.895812

In your code:

locations$Lon <- as.numeric(as.character(locations$Lon))
locations$Lat <- as.numeric(as.character(locations$Lat))