1

I wish to extract raster values based on a list of coordinates. I’ve found online some scripts that include coordinates(), SpatialPoints(), crs() and spTransform() and other that don’t. Could someone kindly explain if script 1 or script 2 is correct and why? Thank you very much!

SCRIPT 1

sites <- read.csv("df.csv")
coordinates(sites)= ~ Longitude+ Latitude 
mypoints = SpatialPoints(sites,proj4string = CRS("+init=epsg:4326"))
myproj = CRS(myraster)
points.proj = spTransform(mypoints, myproj)
myvalues = extract(myraster, points.proj)

SCRIPT 2

sites <- read.csv("df.csv")
myvalues = extract(myraster, cbind(sites$Longitude, y=sites$Latitude), df=TRUE, method='simple', cellnumbers=T)
2
  • can you please paste in the results of head(sites) so we can see the types and structure of the .csv data you're starting with?
    – dbo
    Jul 26, 2019 at 6:42
  • Next time please include some example data like I do in my answer Jul 26, 2019 at 15:15

1 Answer 1

3

Either could be correct. With RasterLayer r and data.frame sites you can do

v <- extract(r, sites[, c("Longitude", "Latitude")])

Under the assumption that "Longitude" and "Latitude" are variables in sites.

However that only works when r also has a ("Longitude", "Latitude") coordinate reference system. That may not be the case. Consider this RasterLayer

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r

#class      : RasterLayer 
#dimensions : 115, 80, 9200  (nrow, ncol, ncell)
#resolution : 40, 40  (x, y)
#extent     : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#crs        : +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +no_defs 
#source     : C:/soft/R/R-3.6.1/library/raster/external/test.grd 
#names      : test 
#values     : 128.434, 1805.78  (min, max)

The crs is "sterea ..." and the extent "178400, 181600, ...) shows that the coordinates are clearly not longitude and latitude (they are expressed in meters away from the origin of the crs.)

In this case, you might have a point in the area covered by r

site <- data.frame(Longitude=5.745039, Latitude=50.96254)

But extract returns NA because the crs do not match

extract(r, site)
#     [,1]
#[1,]   NA

So we do

pts <- SpatialPoints(site)
crs(pts) <- "+proj=longlat +datum=WGS84"
rcrs <- crs(r)
ptrans <- spTransform(pts, rcrs)

And now it works

extract(r, ptrans)
#1529.66 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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