[R-sig-Geo] [raster newbie] projectRaster hangs regridding global latlong -> CONUS LCC
Tom Roche
Tom_Roche at pobox.com
Wed Oct 3 17:34:58 CEST 2012
summary: Robert Hijmans' advice fixed my hang, but now I'm getting
very wrong output--i.e., input variables are not preserved. How to
get the appropriate output schema (e.g., data variables)?
details:
https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016215.html
>> My current [NetCDF] 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
condensed from original (see link)
>> > dimensions:
>> > lon = 360 ;
>> > lat = 180 ;
>> > time = UNLIMITED ; // (1 currently)
>> > variables:
>> > double lon(lon) ;
>> > double lat(lat) ;
>> > double time(time) ;
>> > double emi_n2o(time, lat, lon) ;
>> 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
https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20121002/b9eab9cd/attachment.pl
> If you do not have a template ['to' raster], 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
> # [input] 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')
That fixed the hang problem (see previous links) but the output is
very wrong:
following wrapped for email, and comments added:
> library(raster)
Loading required package: sp
raster 2.0-05 (19-June-2012)
> library(maptools)
Loading required package: foreign
Loading required package: lattice
Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in
maptools depend on gpclib, which has a restricted licence.
It is disabled by default; to enable gpclib, type gpclibPermit()
> gpclibPermit()
[1] TRUE
> library(rgdal)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.6.2, released 2009/07/31
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.6.0, 21 Dec 2007, [PJ_VERSION: 450]
Path to PROJ.4 shared files: (autodetected)
> in.fp <- './GEIA_N2O_oceanic.nc'
> out.fp <- './GEIA_N2O_oceanic_regrid.nc'
# note retry out.crs.2 below
> out.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x=_0=-2556 +y_0=-1728'
# make template raster
data(wrld_simpl)
x <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]
us <- spTransform(x, CRS(out.crs))
out.raster.template <- raster(us)
res(out.raster.template) <- 100000
> in.raster <- raster(in.fp, varname='emi_n2o', band=3)
Loading required package: ncdf
> in.raster[] <- runif(ncell(in.raster))
> 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 +ellps=WGS84 +towgs84=0,0,0
values : in memory
layer name : N2O.emissions
min value : 7.335795e-06
max value : 0.9999727
z-value : 0
# try writing output: fast but very wrong
> system(sprintf("date"))
Wed Oct 3 10:33:08 EDT 2012
> out.raster <-
+ projectRaster(
+ from=in.raster, to=out.raster.template, method='bilinear',
+ overwrite=TRUE, progress='window', format='CDF', filename=out.fp)
> system(sprintf("date"))
Wed Oct 3 10:33:09 EDT 2012
> system(sprintf("ls -alh"))
total 544K
drwxr-xr-x 2 roche mod3 65 Oct 3 10:33 .
drwxr-xr-x 15 roche mod3 4.0K Oct 2 15:25 ..
-rw-r--r-- 1 roche mod3 512K Oct 2 15:28 GEIA_N2O_oceanic.nc
-rw-r--r-- 1 roche mod3 26K Oct 3 10:33 GEIA_N2O_oceanic_regrid.nc
# compare to input netCDF schema
> system(sprintf("ncdump -h %s", out.fp))
netcdf GEIA_N2O_oceanic_regrid {
dimensions:
northing = 93 ;
easting = 66 ;
variables:
double northing(northing) ;
northing:units = "meter" ;
double easting(easting) ;
easting:units = "meter" ;
float variable(easting, northing) ;
variable:units = ;
variable:missing_value = -3.4e+38f ;
variable:_FillValue = -3.4e+38f ;
variable:long_name = "variable" ;
variable:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +y_0=-1728 +ellps=WGS84" ;
variable:projection_format = "PROJ.4" ;
variable:min = -0.01560439f ;
variable:max = 0.9809174f ;
// global attributes:
:Conventions = "CF-1.4" ;
:created_by = "R, packages ncdf and raster (version 2.0-05)" ;
:date = "2012-10-03 10:33:09" ;
}
# retry with different PROJ.4 string, omitting eastings and northings
system(sprintf("rm %s", out.fp))
# compare to out.crs above
out.crs.2 <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97'
# create new template raster
us.2 <- spTransform(x, CRS(out.crs.2))
out.raster.template.2 <- raster(us.2)
res(out.raster.template.2) <- 100000
# create new output: no change in output schema :-(
> system(sprintf("ls -alh"))
total 516K
drwxr-xr-x 2 roche mod3 32 Oct 3 10:42 .
drwxr-xr-x 15 roche mod3 4.0K Oct 2 15:25 ..
-rw-r--r-- 1 roche mod3 512K Oct 2 15:28 GEIA_N2O_oceanic.nc
> system(sprintf("date"))
Wed Oct 3 10:43:29 EDT 2012
> out.raster.2 <-
+ projectRaster(
+ from=in.raster, to=out.raster.template.2, method='bilinear',
+ overwrite=TRUE, progress='window', format='CDF', filename=out.fp)
> system(sprintf("date"))
Wed Oct 3 10:43:30 EDT 2012
> system(sprintf("ls -alh"))
total 544K
drwxr-xr-x 2 roche mod3 65 Oct 3 10:43 .
drwxr-xr-x 15 roche mod3 4.0K Oct 2 15:25 ..
-rw-r--r-- 1 roche mod3 512K Oct 2 15:28 GEIA_N2O_oceanic.nc
-rw-r--r-- 1 roche mod3 26K Oct 3 10:43 GEIA_N2O_oceanic_regrid.nc
> system(sprintf("ncdump -h %s", out.fp))
netcdf GEIA_N2O_oceanic_regrid {
dimensions:
northing = 93 ;
easting = 66 ;
variables:
double northing(northing) ;
northing:units = "meter" ;
double easting(easting) ;
easting:units = "meter" ;
float variable(easting, northing) ;
variable:units = ;
variable:missing_value = -3.4e+38f ;
variable:_FillValue = -3.4e+38f ;
variable:long_name = "variable" ;
variable:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +ellps=WGS84" ;
variable:projection_format = "PROJ.4" ;
variable:min = -0.01560439f ;
variable:max = 0.9809174f ;
// global attributes:
:Conventions = "CF-1.4" ;
:created_by = "R, packages ncdf and raster (version 2.0-05)" ;
:date = "2012-10-03 10:43:29" ;
}
By contrast, I expect output schema appropriate for the input, merely
regridded for the output domain
https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
e.g.,
dimensions:
lon = 459 ; // I should probably rename, e.g. 'COL'
lat = 299 ; // I should probably rename, e.g. 'ROW'
time = UNLIMITED ; // (1 currently)
variables:
double lon(lon) ;
double lat(lat) ;
double time(time) ;
double emi_n2o(time, lat, lon) ;
How to tell projectRaster() what I want? I'm guessing I need to create
a better output template (i.e. argument='to'); how to do that?
your assistance is appreciated, Tom Roche <Tom_Roche at pobox.com>
More information about the R-sig-Geo
mailing list