[GIS] projectRaster() brick from NetCDF changes all values to NA

coordinate systemnetcdfrrasterrgdal

I have loaded a NetCDF brick like so…

library(raster)
library(rgdal)
library(ncdf4)
nctemp <- brick("C:/nco/out.nc")

Which has the following dimensions…

class       : RasterBrick 
dimensions  : 103, 200, 20600, 1680  (nrow, ncol, ncell, nlayers)
resolution  : 0.125, 0.125  (x, y)
extent      : 244, 269, 40, 52.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\nco\out.nc 
names       : X1960.01.16, X1960.02.14, X1960.03.15, X1960.04.14, 
X1960.05.15, X1960.06.14, X1960.07.15, X1960.08.15, X1960.09.14, 
X1960.10.15, X1960.11.14, X1960.12.15, X1961.01.16, X1961.02.14, 
X1961.03.16, ... 
Date        : 1960-01-16, 2099-12-16 (min, max)
varname     : tmin

Further, print(nctemp) returns the following…

2 variables (excluding dimension variables):
double time_bnds[bnds,time]   
double tmin[lon,lat,time]   
    long_name: Downscaled Minimum Temperature in Degrees Celsius
    units: degreesC
    _FillValue: 1e+30
    missing_value: 1e+30

4 dimensions:
lon  Size:200
    standard_name: longitude
    long_name: longitude
    units: degrees_east
    axis: X
lat  Size:103
    standard_name: latitude
    long_name: latitude
    units: degrees_north
    axis: Y

Ultimately, I have a shapefile that I want to use to summarize the pixels. I want to project the raster brick to Albers Equal Area, like so…

nc <- projectRaster(nctemp,crs="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

But when I attempt this or any other transformation, all of the values are changed to NA. I tried to follow along with the suggestions in this question:
https://stackoverflow.com/questions/38356326/projectraster-returns-all-na-values

…and moreover, the linked question on GIS Stack Exchange. However, when I project my shapefile to Albers and then change the extent of the raster to match my shapefile, the projection fails and gives me the error of missing values. Some other possible clues to something weird going on…

  • it can't rotate (raster::rotate gives me an error saying it doesn't "look like an appropriate object for this function")
    a function like raster(inputfile, varname="lat")doesn't work; when I try that,
  • it says that "lat" is not a variable (though it's a dimension for the actual variable of interest, "tmin" in this case, which seems to work for other intents and purposes).

What am I missing?

Best Answer

Your X coordinates need to be -180 to 180. I can illustrate this with a single layer raster:

This has the same extent and res as your brick:

r = raster(xmn=244,xmx=269,ymn=40,ymx=52.875, res=c(0.125,0.125))
r[] = runif(ncell(r))

Now project to your desired CRS:

rp = projectRaster(r,crs="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

> rp
class       : RasterLayer 
dimensions  : 124, 233, 28892  (nrow, ncol, ncell)
resolution  : 9600, 13800  (x, y)
extent      : -1747187, 489613, 1813604, 3524804  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : NA, NA  (min, max)

And the values are NA...

Create it with X coordinates minus 360:

> r = raster(xmn=244-360,xmx=269-360,ymn=40,ymx=52.875, res=c(0.125,0.125))
> r[] = runif(ncell(r))

Project:

> rp = projectRaster(r,crs="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
> rp
class       : RasterLayer 
dimensions  : 124, 233, 28892  (nrow, ncol, ncell)
resolution  : 9600, 13800  (x, y)
extent      : -1747187, 489613, 1813604, 3524804  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : -0.1737072, 1.299517  (min, max)

And now we have real values in rp. Plotting it looks good too.

Note you can reassign an extent to a raster without modifying the data:

> extent(r)
class       : Extent 
xmin        : 244 
xmax        : 269 
ymin        : 40 
ymax        : 52.875 
> extent(r)=c(xmin(r)-360, xmax(r)-360, ymin(r), ymax(r))
> extent(r)
class       : Extent 
xmin        : -116 
xmax        : -91 
ymin        : 40 
ymax        : 52.875
Related Question