在R中将LAT/LON坐标绘制在geotiff上

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()函数从特定位置获取栅格文件中的值。

然而，我的纬度/经度点积累在地图的中心。我错在哪里？非常感谢任何帮助或意见！

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)

- Larusson

---

你的tiff文件不是地理投影，你必须先重新投影它。 - Geo-sp

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"。之后，您可以使用spTransform将数据转换为您想要的crs。

- Robert Hijmans

---

太好了！这个可行...感谢详细的解释。现在我只是在绘图方面有些困难，如果我尝试绘制 plot(r)，我会得到以下错误信息 Error in colnames<-(*tmp*, value = "asi.AMSR2.n6250.20150101.v5") : length of 'dimnames' [2] not equal to array extent。 - Larusson

Answer 2

你可以使用projectRaster来重新投影光栅文件。然后，你可以将点叠加在上面并提取值。

newproj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
p <- projectRaster(r1, newproj)

- Geo-sp