[R-sig-Geo] how to plot Raster*? problems with raster::plot, fields::image.plot
Tom Roche
Tom_Roche at pobox.com
Thu Oct 25 01:30:08 CEST 2012
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 at pobox.com>
More information about the R-sig-Geo
mailing list