根据坐标列表提取栅格值 - spTransform？

Question

我希望根据坐标列表提取栅格值。我在网上找到了一些包含 coordinates()、SpatialPoints()、crs() 和 spTransform() 以及其他不包含的脚本。有人可以解释脚本 1 或脚本 2 是否正确，为什么？非常感谢！

脚本 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)

脚本 2

sites <- read.csv("df.csv")
myvalues = extract(myraster, cbind(sites$Longitude, y=sites$Latitude), df=TRUE, method='simple', cellnumbers=T)

Answer 1

两者都可能是正确的。使用 RasterLayer r 和 data.frame sites 你可以做

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

假设"Longitude"和"Latitude"是sites中的变量。

然而，只有当 r 也有一个 ("Longitude", "Latitude") 坐标参考系统时才有效。情况可能并非如此。考虑这个 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)

crs是"sterea ..."，extent "178400, 181600, ...) 说明坐标明显不是经纬度（以距离原点的米数表示） crs.)

在这种情况下，您可能在 r

覆盖的区域中有一个点

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

但是extractreturnsNA因为crs不匹配

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

我们也是

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

现在可以使用了

extract(r, ptrans)
#1529.66

根据坐标列表提取栅格值 - spTransform？

Extract raster value based on list of coordinates - spTransform?

r

raster

coordinates

r-raster

sp