我有一张海冰图作为“geotiff”。最终目标是提取特定纬度/经度坐标处的海冰浓度。
geotiff 可以在以下位置找到:http://www.iup.uni-bremen.de:8084/amsr2data/asi_daygrid_swath/n6250/2015/jan/asi-AMSR2-n6250-20150101-v5.tif
我想要做的是使用 raster() 加载 geotiff,然后用我的位置覆盖它,然后使用 extract() 函数从特定位置的光栅文件中获取值。
然而,我的纬度/经度点累积在 map 的中心。我哪里出错了?非常感谢任何帮助或输入!
library(raster)
library(sp)
r1 = raster("test.tif")
##check plot
plot(r1)
## check projection
projection(r1)
mydf <- structure(list(longitude = rep(22,7), latitude = seq(60,90,5)),.Names = c("longitude","latitude"), class = "data.frame", row.names = c(NA, -7L))
### Get long and lat from data.frame.
xy <- mydf[,c(1,2)]
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
proj4string = CRS("+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
points(spdf, col="red", lwd=2)
##extract RGB values as reference for sea-ice concentration
seaice_conc = extract(r1, spdf)
最佳答案
Geo-sp 的解决方案可行,但不是最佳的(缓慢且不精确)。您应该始终(重新)投影您的矢量(在这种情况下为点)数据而不是您的栅格数据。投影栅格数据会更改值,而矢量数据则不然。投影栅格数据的计算量也更大。
因此,你应该做这样的事情:
library(raster)
r <- raster("asi-AMSR2-n6250-20150101-v5.tif")
crs(r)
# CRS arguments:
# +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
df <- data.frame(longitude = rep(22,7), latitude = seq(60,90,5), ID=1:7)
spdf <- SpatialPointsDataFrame(coords = df[,1:2], data = df,
proj4string = CRS("+proj=longlat +datum=WGS84"))
library(rgdal)
p <- spTransform(spdf, crs(r))
extract(r, p)
需要注意的是,您犯了一个非常常见的错误:
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf, proj4string =
CRS("+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m"))
您创建了一个
SpatialPointsDataFrame
并分配了错误的坐标引用系统 (crs)。您必须分配实际与您的数据匹配的 crs,即 "+proj=longlat +datum=WGS84"
。之后,您可以将数据转换为您想要的 crs(使用 spTransform
)。关于r - 在 R 中的 geotiff 上绘制 LAT/LON 坐标,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/36453084/