CODEWITHSUNDEEP
rcoordinatesrasterr-rasterr-sp
Cecile
Cecile

Reputation: 535

Extract raster value based on list of coordinates - spTransform?

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)

Upvotes: 1

Views: 482

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47146

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

Upvotes: 3

Related Questions