Tom,

If you do not provide a "to" argument, projectRaster guesses one, and this
is not entirely fool-proof. I think that in this case it created a big
task, and that is simply taking a very long time. Adding an argument like "
progress='window' " can be helpful. projectRaster is 4 times faster om my
PC with 8 cores when I first run beginCluster().

More reasonable would be to set up a reasonable 'to' raster. If you do not
have a template, you can create one like this:

library(raster)
library(maptools)
library(rgdal)

data(wrld_simpl)
x <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]
out.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556
+y_0=-1728'

us <- spTransform(x, CRS(out.crs))
r <- raster(us)
res(r) <- 100000
# or first further reduce to your area of interest by:
# plot(us)
# e <- drawExtent() # click on map twice to draw box
# r <- raster(e, crs=out.crs)
# res(r) <- 100000

# data similar to yours
in.raster <- raster()
in.raster[] <- runif(ncell(in.raster))

out <- projectRaster(from=in.raster, to=r, method='bilinear',
overwrite=TRUE, progress='window')
plot(out)
plot(us, add=TRUE)


Best, Robert


On Tue, Oct 2, 2012 at 2:29 PM, Tom Roche <Tom_Roche@pobox.com> wrote:

>
> 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@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@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@pobox.com>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

