Extract values from raster from list of text file of coordinates in r

Question

I have several raster files that i am trying to extract values from that come from some lat longs found in a text file. I've been successful with raster::extract in the past but seem to be running into a new problem here. The extract function returns NAs despite the projections seeming to line up.

Here is an (almost) reproducible example that shows my issue.

#reproducable example for SO
#
library(leaflet)
library(rgeos)
library(rgdal)
library(raster)
testRast <- raster("FranceRasters/FRA_STU_EU_S_CLAY.rst")
# crs string for ETRS89
crs(testRast) <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

testRast
# class       : RasterLayer 
# dimensions  : 1109, 1076, 1193284  (nrow, ncol, ncell)
# resolution  : 1000, 1000  (x, y)
# extent      : 3209000, 4285000, 2026000, 3135000  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# data source : D:\GIS\STU_EU_Layers\FranceRasters\FRA_STU_EU_S_CLAY.rst 
# names       : FRA_STU_EU_S_CLAY 
# values      : 0, 61  (min, max)

#Small sample of lat ongs that need sampling
Coordinate.sample <-  structure(list(Longitude = c(7.62648877, 7.63023291, 7.37533534, 
                                                   7.37484746, 7.59758433, 7.46470834, 7.46133708, 7.46632616, 7.45291197, 
                                                   7.52007514), Latitude = c(48.56788166, 48.59294834, 48.65367653, 
                                                                             48.65288703, 48.78729009, 48.71844903, 48.71878129, 48.71828592, 
                                                                             48.71898041, 48.52164554)), .Names = c("Longitude", "Latitude"
                                                                             ), row.names = c(NA, 10L), class = "data.frame")

#convert to spatial points
Coordinate.sample <- SpatialPoints(Coordinate.sample)

projection(Coordinate.sample) <- projection(testRast)
#show projections line up
leaflet() %>% addTiles() %>% addRasterImage(testRast) %>% addCircleMarkers(data=Coordinate.sample)

#extraction leads to NAs
raster::extract(testRast,Coordinate.sample)
      [,1]
 [1,]   NA
 [2,]   NA
 [3,]   NA
 [4,]   NA
 [5,]   NA
 [6,]   NA
 [7,]   NA
 [8,]   NA
 [9,]   NA
[10,]   NA

As you can see the extraction has failed. Does anyone have any ideas why it has failed?

Here is an image of the leaflet() output to show that the coords and raster are in the right projection.

Answer 1

You are overwriting the projection of Coordinate.sample instead of transforming it.

First set the original coordinate system when creating SpatialPoints, in this case it seems to be WGS84:

Coordinate.sample <- SpatialPoints(Coordinate.sample, proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

and then change the line

projection(Coordinate.sample) <- projection(testRast)

to

Coordinate.sample <- spTransform(Coordinate.sample, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))

and it should work. This will transform your coordinates to the specified projection.