1

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?

  • Another curious thing that makes me wonder what's up with my raster is that I can't rotate it... – setophaga Jun 30 '17 at 17:31
2

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 
2

I would not project this raster, but use

extract(shift(nctemp, x = -360), shape)

where "shape" is the Spatial*DataFrame read in from shapefile. Whether this will work depends on the correctness of the object's CRS metadata, but to be sure try

plot(shift(nctemp[[1]], x = -360), addfun = function() plot(spTransform(shape, projection(nctemp)), add = TRUE))

You can use

nctemp <- shift(nctemp, x = -360)

when you first create it or after, as long as you only do it once, at the right moment. :)

Extracting with a reprojected shapefile is far higher fidelity than warping a huge raster, though it does depend ultimately on the topological layout of your shapes and the nature of your workflow.

  • Agree with your reasoning that it's better to project the shape file than warp the raster; I think shift only takes non-negative integers though. So will combine your wisdom with the answer above to just manipulate the extent of the raster and project the shape file to match the raster's projection! rdocumentation.org/packages/data.table/versions/1.10.4/topics/… – setophaga Jul 5 '17 at 15:59
  • 1
    raster::shift will definitely take negative integers, if it doesn't let me know and we'll force it into submission one way or another – mdsumner Jul 6 '17 at 5:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.