r - 在 R 中的 geotiff 上绘制 LAT/LON 坐标

Question

我有一张海冰图作为“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)

Answer 1

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)。