如何在同一投影上获取此栅格和此 shapefile？

Question

我有一个在 NAD83（硬）/CA Albers 中投影的 10 个 CA 县的 shapefile。我在 WGS84/WGS84 中有整个美国的光栅（温度的 netCDF 文件）。我想使用 shapefile 来裁剪栅格。我知道我需要先让他们在同一个 datum/projection 上。但我尝试使用 raster::projectRaster() 重新投影栅格。那失败了（因为数据消失了）。因此，我尝试使用 sp::spTransform() 重新投影 shapefile。这也失败了（因为数据不重叠）。我搜索了 Whosebug 但没有看到任何似乎有帮助的东西。我没有收到错误，但是 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)

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

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

如何在同一投影上获取此栅格和此 shapefile？

How do I get this raster and this shapefile on the same projection?

r

projection

raster

r-raster