Как мне получить этот растр и этот шейп-файл в одной проекции?

Question

У меня есть шейп-файл из 10 округов CA, спроецированный в NAD83 (Hard) / CA Albers.У меня есть растр (файл температуры netCDF) для всего США, спроецированного в WGS84 / WGS84.Я хочу использовать шейп-файл для обрезки растра.Я знаю, что сначала мне нужно получить их в одном и том же датуме / проекции.Но я попытался перепроектировать растр с помощью raster :: projectRaster ().Это не удалось (как в данных исчезли).Тогда я попытался перепроектировать шейп-файл вместо этого, используя sp :: spTransform ().Это также не удалось (так как данные не перекрываются).Я искал через stackoverflow, но не увидел ничего, что могло бы помочь.Я не получаю сообщение об ошибке, но projectRaster не работает, и перепроектирование шейп-файла с использованием spTransform не дает желаемого результата.Я чувствую, что здесь происходит что-то конкретное, например, преобразование из WGS84 в NAD83 или загрузка растра с использованием raster() - проблема ... но опять же, это может быть что-то глупое, что я упускаю!=)

мой шейп-файл и растр здесь: 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)

Обычно лучше преобразовывать векторные, а не растровые данные, если вы можете избежать этого.