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

【问题标题】：Extract raster value based on list of coordinates - spTransform?根据坐标列表提取栅格值 - spTransform？
【发布时间】：2019-12-04 10:50:39
【问题描述】：

我希望根据坐标列表提取栅格值。我在网上找到了一些脚本，包括坐标()、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)

【问题讨论】：

能否请您粘贴head(sites) 的结果，以便我们查看您开始使用的.csv 数据的类型和结构？
下次请像我在回答中那样包含一些示例数据

标签： r coordinates raster r-raster sp

【解决方案1】：

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

v <- extract(r, sites[, c("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 ...”，范围“178400, 181600, ...) 表明坐标显然不是经度和纬度（它们以距离 crs 原点的米数表示。 )

在这种情况下，您可能在r 覆盖的区域内有一个点

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

但是 extract 返回 NA 因为 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

【讨论】：