Tom, I do not know why you say the output is wrong. Why you would expect that input variables would have to be "preserved" (copied to the output file)? For 'raster' there is input data and output data, the input and output files (formats, attributes) are not related in any way.
You could use the ncdf(4) package to further manipulate the data and move them from one file to another. 'raster' could also conceivably offer more functionality to specify the properties of the ncdf file you want (see the additional arguments for writing to a ncdf file in ?writeRaster), and suggestions are always welcome. Robert On Wed, Oct 3, 2012 at 8:34 AM, Tom Roche <tom_ro...@pobox.com> wrote: > > 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_ro...@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]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
