[R-sig-Geo] [raster newbie] projectRaster hangs regridding global latlong -> CONUS LCC

[Tom Roche]

summary: I'm debugging an apparent hang of raster::projectRaster when
projecting the values from an input global lon-lat-gridded emissions
inventory (EI) to an output over the contiguous US (CONUS) gridded
Lambert conformal conic (LCC). I'm new to 'raster' and relatively new
to geomatics and geospatial transforms. I suspect the problem involves
failure to provide either argument='to' or argument='res', but that
raises other questions.

details:

I'm learning package=raster (more questions to follow!) by RTFM,
notably

http://cran.r-project.org/web/packages/raster/raster.pdf
http://cran.r-project.org/web/packages/raster/vignettes/Raster.pdf

and the CLI help, due to a truncation problem

https://stat.ethz.ch/pipermail/r-sig-geo/2012-September/016179.html

I'm working with netCDF. My current input

https://github.com/downloads/TomRoche/ioapi-hack-R/GEIA_N2O_oceanic.nc

is gridded lon-lat @ 1°x1° with global extents:

$ ncdump -h ./GEIA_N2O_oceanic.nc 
> netcdf GEIA_N2O_oceanic {
> dimensions:
>       lon = 360 ;
>       lat = 180 ;
>       time = UNLIMITED ; // (1 currently)
> variables:
>       double lon(lon) ;
>               lon:units = "degrees_east" ;
>               lon:long_name = "longitude" ;
>               lon:comment = "center_of_cell" ;
>       double lat(lat) ;
>               lat:units = "degrees_north" ;
>               lat:long_name = "latitude" ;
>               lat:comment = "center_of_cell" ;
>       double time(time) ;
>               time:units = "year" ;
>               time:long_name = "time" ;
>               time:calendar = "none: emissions attributed to no particular year" ;
>       double emi_n2o(time, lat, lon) ;
>               emi_n2o:units = "ton N2O-N/yr" ;
>               emi_n2o:missing_value = -999. ;
>               emi_n2o:long_name = "N2O emissions" ;
>               emi_n2o:_FillValue = -999.f ;

In order to combine this with other inputs, I need to regrid this to
the "AQMEII North American" domain

https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA 

for which I believe the PROJ.4 string is

+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556 +y_0=-1728

(correct? I'm pretty sure about everything but the easting/northings,
which I'm not sure if I even need.) Note the extents of AQMEII-NA are
over CONUS, not global.

So I tried to do the following, on a cluster on which I don't have root
(i.e., I can't just upgrade stuff myself):

me at it:/project/inf14w/me/regrid_raster $ R
R version 2.14.1 (2011-12-22)
...
> library(raster)
Loading required package: sp
raster 2.0-05 (19-June-2012)
> library(ncdf4)
> in.fp <- './GEIA_N2O_oceanic.nc'
> system(sprintf("ls -alh %s", in.fp))
-rw-r--r-- 1 me mod3 512K Oct  2 15:28 ./GEIA_N2O_oceanic.nc
> in.raster <- raster(in.fp, varname='emi_n2o', band=3)
Loading required package: ncdf
> in.raster
class       : RasterLayer 
dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
values      : /project/inf14w/me/regrid_raster/GEIA_N2O_oceanic.nc 
layer name  : N2O.emissions 
z-value     : 0 
zvar        : emi_n2o 

> out.fp <- './GEIA_N2O_oceanic_regrid.nc'
> out.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556 +y_0=-1728'
> out.raster <- 
+   projectRaster(
+     from=in.raster, crs=out.crs, method='bilinear', format='CDF', filename=out.fp)

This ran for a few minutes (with no stdout from projectRaster(...)),
at which point I observed something like

me at it:/project/inf14w/me/regrid_raster $ ls -alh
total 147M
drwxr-xr-x  2 me mod3   65 Oct  2 16:02 .
drwxr-xr-x 15 me mod3 4.0K Oct  2 15:25 ..
-rw-r--r--  1 me mod3 512K Oct  2 15:28 GEIA_N2O_oceanic.nc
-rw-r--r--  1 me mod3 146M Oct  2 16:08 GEIA_N2O_oceanic_regrid.nc

I let it run a few minutes more before killing (with ^C in R); folder
contents were unchanged. So I'm wondering, what would cause this hang?
and how to debug?

I tend to suspect either

1 the absence of argument=to, created by a call to 

projectExtent(object, crs)

2 the absence of argument=res, given the absence of argument=to

per the doc

?projectRaster
> Usage:

>   projectRaster(from, to, res, crs, method="bilinear", filename="", ...) 
>   projectExtent(object, crs)
     
> Arguments:

>   from: Raster* object

>     to: Raster* object with the parameters to which 'from' should be
>         projected

>    res: Single or (vector of) two numerics. To, optionally, set the
>         output resolution if 'to' is missing

I have files with the form of the outputs I need (i.e., netCDF gridded
for AQMEII-NA), but they all use IOAPI

http://www.baronams.com/products/ioapi/

(a netCDF wrapper format) which could be confusing to projectRaster.
I'll try that later, in the absence of alternatives, since I don't
know what to give as 'res', either: I know from

https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA

the output will have ncol=459 and nrow=299, but I don't know the
vertical and horizontal extents of the output (which is not global) in
degrees (which I'm guessing is the numerator of the resolution, with the
denominator being the counts above), and the domain is not rectangular
in lon-lat (see the image just above the previous link). Am I
misunderstanding how to define 'res'?

And what I should do to debug the hang?

Your assistance is appreciated, Tom Roche <Tom_Roche at pobox.com>