[R-sig-Geo] [raster] unexpected write behavior with writeRaster, projectRaster

[Tom Roche]

Noticing an unwelcome similarity between writeRaster and projectRaster:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016523.html
> Here I show the raster object, I then define the name, and then I
> write it to disk. You can see that the layer name switches back to
> "layer".

> > comb
> class       : RasterLayer
> dimensions  : 7200, 7200, 51840000  (nrow, ncol, ncell)
> resolution  : 0.008333333, 0.008333333  (x, y)
> extent      : -150, -90, 30, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source :
> /private/var/folders/4t/f1r24pf90d724kpr9m7bpl640000gp/T/R_raster_tmp/Pascal/raster_tmp_2012-10-26_101222_71239.grd

> names       : layer
> values      : -290, 239  (min, max)

> > names(comb)<-'bio1'
> > comb
> class       : RasterLayer
> dimensions  : 7200, 7200, 51840000  (nrow, ncol, ncell)
> resolution  : 0.008333333, 0.008333333  (x, y)
> extent      : -150, -90, 30, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source :
> /private/var/folders/4t/f1r24pf90d724kpr9m7bpl640000gp/T/R_raster_tmp/Pascal/raster_tmp_2012-10-26_101222_71239.grd

> names       : bio1
> values      : -290, 239  (min, max)

> > writeRaster(comb,filename='/Users/Pascal/Desktop/bio1.grd')
> class       : RasterLayer
> dimensions  : 7200, 7200, 51840000  (nrow, ncol, ncell)
> resolution  : 0.008333333, 0.008333333  (x, y)
> extent      : -150, -90, 30, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : /Users/Pascal/Desktop/bio1.grd
> names       : layer
> values      : -290, 239  (min, max)

I see a similar bug with projectRaster: the output CRS is longlat
(presumably the default), even if *both* the extents object ('to')
used to create it is projected (i.e., has a non-longlat CRS), *and*
I pass the desired CRS to projectRaster(...), although the crs slot in
the written netCDF is correct:

> out.crs
[1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000"
> template.in.raster <- raster(template.in.fp, varname=data.var.name, band=template.band)
> template.raster <- projectExtent(template.in.raster, crs=out.crs)
Warning message:
In projectExtent(template.in.raster, crs = out.crs) :
  158 projected point(s) not finite

Is that "projected point(s) not finite" warning important?
Probably not, per Hijmans.

> template.raster
class       : RasterLayer 
dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
resolution  : 53369.55, 56883.69  (x, y)
extent      : -14802449, 9694173, -6258782, 10749443  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000 +ellps=WGS84 

That is the correct CRS. Just to be sure,

> template.raster at crs
CRS arguments:
 +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000
+y_0=-1728000 +ellps=WGS84 

> out.raster <-
+   projectRaster(
+     from=in.raster, to=template.raster, method='bilinear', crs=out.crs,
+     overwrite=TRUE, progress='window', format='CDF',
+     # args from writeRaster
+     NAflag=-999.0,  # match emi_n2o:missing_value,_FillValue (TODO: copy)
+     varname=data.var.name, 
+     varunit='ton N2O-N/yr',
+     longname='N2O emissions',
+     xname='COL',
+     yname='ROW',
+     filename=out.fp)
> out.raster
class       : RasterLayer 
dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
resolution  : 53369.55, 56883.69  (x, y)
extent      : -14802449, 9694173, -6258782, 10749443  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /path/to/GEIA_N2O_oceanic_regrid.nc 
names       : N2O.emissions 
zvar        : emi_n2o 

That CRS is incorrect, but ...

> system(sprintf('ncdump -h %s', out.fp))
...
    float emi_n2o(ROW, COL) ;
      emi_n2o:units = "ton N2O-N/yr" ;
      emi_n2o:_FillValue = -999. ;
      emi_n2o:missing_value = -999. ;
      emi_n2o:long_name = "N2O emissions" ;
      emi_n2o:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000 +ellps=WGS84" ;

... that CRS is correct. Inspecting the out.raster CRS,

> out.raster at crs
CRS arguments:
 +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
> out.raster at crs <- CRS(out.crs)
> out.raster
class       : RasterLayer 
dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
resolution  : 53369.55, 56883.69  (x, y)
extent      : -14802449, 9694173, -6258782, 10749443  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000 +ellps=WGS84 
data source : /path/to/GEIA_N2O_oceanic_regrid.nc 
names       : N2O.emissions 
zvar        : emi_n2o 

Unfortunately fixing out.raster's CRS does not fix the problem with
raster::plot-ing it: plot output still has incorrect (global) extents,
as described @

https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016504.html

FWIW, Tom Roche <Tom_Roche at pobox.com>