1

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.

enter image description here

3

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.

4
  • Hi Janina, running your suggestion gives me the error No transformation possible from NA reference system – Michael Bird Oct 16 '17 at 13:04
  • @MichaelBird Oh, sorry, I didn't realise that the coordinate system wasn't set in the first place. Give me a minute and I'll edit my answer. – Janina Oct 16 '17 at 13:06
  • @MichaelBird Edited. – Janina Oct 16 '17 at 13:12
  • Wonderful, that's working perfectly now. Thank you very much! – Michael Bird Oct 16 '17 at 13:13

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.