CODEWITHSUNDEEP
rrasterprojectionr-raster
ESELIA
ESELIA

Reputation: 162

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

I have a shapefile of 10 CA counties projected in NAD83(Hard)/CA Albers. I have a raster (a netCDF file of temperature) for the entire US projected in WGS84/WGS84. I want to use the shapefile to clip the raster. I know that I need to get them on the same datum/projection first. But I've tried re-projecting the raster using raster::projectRaster(). That failed (as in the data disappeared). So then I tried re-projecting the shapefile instead using sp::spTransform(). This also failed (as in the data don't overlap). I've searched through stackoverflow but didn't see anything that seemed to help. I'm not getting an error, but projectRaster is not working and re-projecting the shapefile using spTransform doesn't produce the desired outcome. I feel like there is something specific going on here, like the transformation from WGS84 to NAD83 or loading the raster in using raster() is the problem ... but then again, it could easily be something stupid that I'm missing! =)

my shapefile and raster are here: https://www.dropbox.com/sh/l3b2syzcioeqmyy/AAA5CstBZty4ofOcVFkAumNYa?dl=0

here is my code:

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)

Upvotes: 0

Views: 733

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47071

You can do

library(raster)
library(rgdal)

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

Now look at 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

The horizontal extent is between 234 and 294 degrees. That points to a system with longitudes that start in Greenwich at 0 and continues to 360 (again in Greenwich). Climatologist do that. To go to the more conventional -180 to 180 degrees system:

r <- shift(b, -360)

(if your data had a global extent, you would use raster::rotate instead)

Now, transform the counties to lonlat and show that they overlap.

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

It is generally best to transform vector data, not raster data if you can avoid it.

Upvotes: 2

Related Questions