Tom,

Please provide a clear question and, whenever possible, a simple, concise,
and self contained program that illustrates the problem. Based on what I
see in your message, here is what I think you do:

library(raster)
# replacing your file with example data
in.raster <- raster()
in.raster[] <- 1:ncell(in.raster)
template.raster <- raster(nrow=299, ncol=459, xmn=-14802449, xmx=9694173,
ymn=-6258782, ymx=10749443, crs="+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')

library(maptools)
data(wrld_simpl)
map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
map.us.proj <- spTransform(map.us.unproj, CRS=projection(out.raster,
asText=F))
plot(out.raster)
plot(map.us.proj, add=TRUE)

# I take it you want a smaller extent. Then use a template with a better
extent
# to get this drew an extent on the map, using
# x <- crop(out.raster, drawExtent())

temp2 <- raster(nrow=70, ncol=117, xmn=-5622887, xmx=621350, ymn=-3642132,
ymx=339726.7, crs="+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97
+x_0=-2556000 +y_0=-1728000 +ellps=WGS84")
out2 <- projectRaster(from=in.raster, to=temp2, method='bilinear')

plot(out2)
plot(map.us.proj, add=TRUE)

# better now? If not, why not?

Note that this script does not write to ncdf (or any other) files at this
point as it is best to solve coding problems one at a time. It might also
be that this works for the example, but not for your file, but we need to
know that before we can move on.

Robert



On Wed, Oct 24, 2012 at 4:30 PM, Tom Roche <Tom_Roche@pobox.com> wrote:

>
> https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016467.html
> >> Could I get one or more spatial-plotting experts to examine defective
> >> output produced when I try to plot a RasterLayer[, as described at]
>
> >>
> http://stackoverflow.com/questions/13005181/rasterplot-problems-unwanted-rotation-extents
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016485.html
> > It would be useful to provide a simple self-contained example (of R
> > code) of what you are trying to achieve.
>
> I have described this completely at the above link, and all the code and
> data is contained in a single git repository. All that's required to run
> it is to clone the repo and run one script.
>
> But I "boil it down" further here, and focus solely on the raster::plot
> problem, on which I have made progress. About the fields::image.plot
> problem I still have no clue; with raster::plot, I fixed the rotation
> problem (which appears to have been due to a former typo in my CRS), but
> am still having the extents problem: the plot of the regridded output
> has much larger extents (global, or close) than does the extents
> RasterLayer used to create it, and therefore to define the regridding.
>
> Here's a simple demonstration:
>
> 1. Open a shell and create a working directory: in bash,
>
> export TEST_DIR='/tmp/projectRasterTest' # or your name here
> rm -fr ${TEST_DIR}   # start fresh
> mkdir -p ${TEST_DIR}
> # sometimes files seem to write to cwd, so move to it just in case
> pushd ${TEST_DIR}
>
> 2. Download netCDF input: in bash,
>
> DATA_INPUT_URI='
> https://github.com/downloads/TomRoche/GEIA_to_netCDF/GEIA_N2O_oceanic.nc'
> DATA_INPUT_FN="$(basename ${DATA_INPUT_URI})"
> export DATA_INPUT_FP="${TEST_DIR}/${DATA_INPUT_FN}"
> if [[ ! -r "${DATA_INPUT_FP}" ]] ; then
>   for CMD in \
>     "wget --no-check-certificate -c -O ${DATA_INPUT_FP} ${DATA_INPUT_URI}"
> \
>   ; do
>     echo -e "$ ${CMD}"
>     eval "${CMD}"
>   done
> fi
> export DATA_VAR_NAME='emi_n2o'  # see ncdump output
> export DATA_INPUT_BAND='3'      # GEIA makes dim=TIME third
>
> 3. Download netCDF used as a template for the extents:
>
> TEMPLATE_INPUT_URI='
> https://github.com/downloads/TomRoche/GEIA_to_netCDF/emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT.nc
> '
> TEMPLATE_INPUT_FN="$(basename ${TEMPLATE_INPUT_URI})"
> export TEMPLATE_INPUT_FP="${TEST_DIR}/${TEMPLATE_INPUT_FN}"
> if [[ ! -r "${TEMPLATE_INPUT_FP}" ]] ; then
>   for CMD in \
>     "wget --no-check-certificate -c -O ${TEMPLATE_INPUT_FP}
> ${TEMPLATE_INPUT_URI}" \
>   ; do
>     echo -e "$ ${CMD}"
>     eval "${CMD}"
>   done
> fi
> export TEMPLATE_INPUT_BAND='1' # `ncks` makes dim=TSTEP first
>
> 4. Setup a few more variables for later use: in bash,
>
> export PDF_VIEWER='xpdf'  # whatever works on your platform
> PDF_FN='projectRasterTest.pdf'
> export PDF_FP="${TEST_DIR}/${PDF_FN}" # path to PDF output
>
> DATA_INPUT_FN_ROOT="${DATA_INPUT_FN%.*}"
> DATA_INPUT_FN_EXT="${DATA_INPUT_FN#*.}"
> DATA_OUTPUT_FN="${DATA_INPUT_FN_ROOT}_regrid.${DATA_INPUT_FN_EXT}"
> export DATA_OUTPUT_FP="${TEST_DIR}/${DATA_OUTPUT_FN}"
>
> 5. From here on, everything is R, so start R in your shell, and load
>    required libraries and data:
>
> library(raster)
> library(ncdf4)
> library(maptools)
> # run if rgeos not available
> # gpclibPermit()
> library(rgdal)
> data(wrld_simpl) # from maptools
>
> 6. Initialize vars, mostly with environment values:
>
> pdf.fp <- Sys.getenv('PDF_FP')
> pdf.er <- Sys.getenv('PDF_VIEWER')
> in.fp <- Sys.getenv('DATA_INPUT_FP')
> in.band <- Sys.getenv('DATA_INPUT_BAND')
> data.var.name <- Sys.getenv('DATA_VAR_NAME')
> template.in.fp <- Sys.getenv('TEMPLATE_INPUT_FP')
> template.band <- Sys.getenv('TEMPLATE_INPUT_BAND')
> out.fp <- Sys.getenv('DATA_OUTPUT_FP')
>
> # This corresponds to the output domain specification @
> #
> https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
> out.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97
> +x_0=-2556000 +y_0=-1728000'
>
> 7. Rasterize the input netCDF:
>
> in.raster <- raster(in.fp, varname=data.var.name, band=in.band)
> 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      : /tmp/projectRasterTest/GEIA_N2O_oceanic.nc
> # layer name  : N2O.emissions
> # z-value     : 0
> # zvar        : emi_n2o
>
> 8. Rasterize the extents netCDF, which should resemble the
>    output domain specification:
>
> template.in.raster <- raster(template.in.fp, varname=data.var.name,
> band=template.band)
> template.raster <- projectExtent(template.in.raster, crs=out.crs)
> 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
>
>    It does!
>
> 9. Regrid the input netCDF using the extents netCDF:
>
> 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
>     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 : /tmp/projectRasterTest/GEIA_N2O_oceanic_regrid.nc
> # names       : N2O.emissions
> # zvar        : emi_n2o
>
>     I'm not sure why the regridded output is still proj=longlat, given
>     that I clearly set a projected CRS on it, and the extents ('to')
>     used to create it has the same CRS.
>
> 10. Examine the regridded netCDF written by projectRaster: looks good,
>     and note that the regridded netCDF knows its CRS (even if its
>     parent RasterLayer does not).
>
> system(sprintf('ncdump -h %s', out.fp))
> # netcdf GEIA_N2O_oceanic_regrid {
> # dimensions:
> #   COL = 459 ;
> #   ROW = 299 ;
> # variables:
> #   double COL(COL) ;
> #       COL:units = "meter" ;
> #       COL:long_name = "COL" ;
> #   double ROW(ROW) ;
> #       ROW:units = "meter" ;
> #       ROW:long_name = "ROW" ;
> #   float emi_n2o(ROW, COL) ;
> #       emi_n2o:units = "ton N2O-N/yr" ;
> #       emi_n2o:missing_value = -999.f ;
> #       emi_n2o:_FillValue = -999.f ;
> #       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" ;
> #       emi_n2o:projection_format = "PROJ.4" ;
> #       emi_n2o:min = 5.9605e-08f ;
> #       emi_n2o:max = 783.92f ;
> #
> # // global attributes:
> #               :Conventions = "CF-1.4" ;
> #               :created_by = "R, packages ncdf and raster (version
> 2.0-12)" ;
>
> typo: should be 'ncdf4'                    ^^^^
>
> 11. Get a map of North America, and project it:
>
> map.us.unproj <-
>   wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
> map.us.proj <- spTransform(map.us.unproj, CRS(out.crs))
>
> 12. Plot the output 3 different ways:
>
> # plot to PDF---------------------------------------------------------
> # repeat for each plot file
> pdf(file=pdf.fp, width=5.5, height=4.25)
>
> # plot page 1: just plain plot----------------------------------------
> plot(out.raster)
> # add the projected CONUS map
> plot(map.us.proj, add=TRUE)
>
> # plot page 2: try setting plot extents-------------------------------
> plot(out.raster, ext=template.raster)
> plot(map.us.proj, add=TRUE)
> # note page 2 is identical to page 1
>
> # plot page 3: try raster::crop---------------------------------------
> template.raster.extent <- extent(template.raster)
> out.raster.crop <-
>   # overwrite the uncropped file
>   crop(out.raster, template.raster.extent, filename=out.fp, overwrite=TRUE)
> plot(out.raster.crop)
> plot(map.us.proj, add=TRUE)
> # note page 3 is identical to page 2 and page 1
>
> # flush the plot device
> dev.off()
>
> 13. Exit R, and examine the plot: in bash,
>
> # After exiting R, show cwd and display output PDF.
> for CMD in \
>   "ls -alht ${TEST_DIR}" \
>   "${PDF_VIEWER} ${PDF_FP} &" \
> ; do
>   echo -e "$ ${CMD}"
>   eval "${CMD}"
> done
>
>     The resultant PDF will have 3 pages, all identical, and all having
>     much greater extents than does
>
> https://github.com/TomRoche/cornbeltN2O/wiki/images/Cmaq_aqmeii_domain.png
>
>     the boundaries of which are the extents of the output domain.
>     How to fix?
>
> 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]]

