如何在同一个投影上获得这个栅格和这个 shapefile？

Posted 2023-02-25

【中文标题】如何在同一个投影上获得这个栅格和这个 shapefile？

【英文标题】：How do I get this raster and this shapefile on the same projection?

【发布时间】：2019-06-15 00:43:59

【问题描述】：

我有一个包含 10 个 CA 县的 shapefile，预计在 NAD83（硬）/CA Albers。我有一个以 WGS84/WGS84 投影的整个美国的栅格（温度的 netCDF 文件）。我想使用 shapefile 来剪辑栅格。我知道我需要先将它们放在相同的基准/投影上。但是我尝试使用 raster::projectRaster() 重新投影光栅。那失败了（因为数据消失了）。因此，我尝试使用 sp::spTransform() 重新投影 shapefile。这也失败了（因为数据不重叠）。我已经通过 *** 进行了搜索，但没有看到任何似乎有帮助的东西。我没有收到错误消息，但 projectRaster 无法正常工作，并且使用 spTransform 重新投影 shapefile 不会产生预期的结果。我觉得这里发生了一些特定的事情，比如从 WGS84 到 NAD83 的转换或使用raster() 加载栅格是问题......但话又说回来，我很容易错过一些愚蠢的事情！ =)

我的 shapefile 和光栅在这里：https://www.dropbox.com/sh/l3b2syzcioeqmyy/AAA5CstBZty4ofOcVFkAumNYa?dl=0

这是我的代码：

library(raster) #for creating rasters from .bil files
library(rgdal) #for reading .bil files and .gdb files
library(ncdf4) #for working with ncdf files
library(sp) #for working with spatial data files

load(my_counties.RData)
myraster <- raster(myraster.nc)
my.crs <- CRS("+init=EPSG:3311") #NAD83(HARN) / California Albers (HARN is high resolution)

newraster <- projectRaster(myraster, res = 6000, crs = my.crs) #raster resolution is 1/16th of a degree

#There is data in the raster.
plot(myraster)

#but none in newraster
plot(newraster)

#Now try re-projecting the shapefile
my.crs2 <- crs(myraster)
newshapefile <- spTransform(my_counties, my.crs2)

#but the data don't overlap
plot(newshapefile); plot(myraster, add = T)

【参考方案1】：

你可以的

library(raster)
library(rgdal)

load("my_counties.RData")
b <- brick("myraster.nc")

现在看b

b
#class      : RasterBrick 
#dimensions : 490, 960, 470400, 365  (nrow, ncol, ncell, nlayers)
#resolution : 0.0625, 0.0625  (x, y)
#extent     : 234, 294, 23.375, 54  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : myraster.nc 
#names      : X2005.01.01, X2005.01.02, X2005.01.03, X2005.01.04, X2005.01.05, X2005.01.06, X2005.01.07, X2005.01.08, X2005.01.09, X2005.01.10, X2005.01.11, X2005.01.12, X2005.01.13, X2005.01.14, X2005.01.15, ... 
#Date       : 2005-01-01, 2005-12-31 (min, max)
#varname    : tasmax 

水平范围在 234 到 294 度之间。这指向一个系统，其经度从格林威治 0 开始，一直到 360 度（同样在格林威治）。气候学家这样做。转到更传统的 -180 到 180 度系统：

r <- shift(b, -360)

（如果您的数据具有全局范围，则可以改用 raster::rotate）

现在，将县转换为 lonlat 并显示它们重叠。

counties <- spTransform(my_counties, crs(r))
plot(r, 1)
lines(counties)

如果可以避免，通常最好转换矢量数据，而不是栅格数据。

【讨论】：

感谢您简洁明了的解释！这也解释了为什么我在尝试 raster::rotate 时遇到错误（“看起来不像这个函数的合适对象”）。在改变范围后，我还能够让 spTransform 处理 shapefile！