[R-sig-Geo] [raster newbie] plot configuration problems
Tom Roche
Tom_Roche at pobox.com
Tue Oct 16 08:00:16 CEST 2012
summary: thanks for your helpful reply
https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016263.html
to my previous post--it seems to have solved the previous problems.
However, when I try to plot the resulting data, I get
http://tinyurl.com/rasterPlotExample_20121016
What am I doing wrong?
details:
Apologies for my delay in responding to
https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016263.html
but untimely commitments intervened. After applying your advice, I now
wish to visualize the data for quality control. I'm attempting to use
raster::plot, and to give it the appropriate extents, but I'm not
getting what I want/expect. Instead,
http://tinyurl.com/rasterPlotExample_20121016
shows the desired view, but very compressed and distorted, plus much
additional material. I'd like to see something with a shape more like
this image
https://github.com/TomRoche/cornbeltN2O/wiki/images/CMAQ_AQMEII_MSFs.png
of the AQMEII-NA domain, which has the same projection (CONUS LCC) and
matrix layout (col=459, row=299) as the data I'm attempting to plot.
So what am I doing wrong?
The following example (quite similar to previous) starts in bash, then
goes to R. It contains some workarounds for another cluster on which
I'm working: please adjust to local requirements.
# begin example-------------------------------------------------------
# do what's required for your R to load ncdf4.so, libnetcdf.so.7
# on my cluster it's
module add netcdf-4.1.2
# and point to whatever displays PDFs there
export PDF_VIEWER='xpdf'
# since I must plot to PDF for visualization, due to X forwarding restrictions
# test folder, made/destroyed
export TEST_DIR='/tmp/projectRasterTest'
rm -fr ${TEST_DIR}
mkdir -p ${TEST_DIR}
# Input data is GEIA marine N2O emissions converted to netCDF:
# see code @ https://github.com/TomRoche/GEIA_to_netCDF
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}"
wget --no-check-certificate -c -O ${DATA_INPUT_FP} ${DATA_INPUT_URI}
export DATA_VAR_NAME='emi_n2o' # see ncdump output
export DATA_INPUT_BAND='3' # GEIA makes dim=TIME third
# output name is variant of input name
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}"
# The template input is a copy of some meteorological data with 52
# "real" datavars. That data is in the IOAPI format, which is in this
# case basically a wrapper around netCDF. (IOAPI can be used with other
# data, and similar wrappers exist for various models.) I removed all
# but one of the datavars (with NCO 'ncks'), and renamed (with NCO
# 'ncrename') that datavar to have the same name as my input datavar
# (and hopefully my output datavar).
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}"
wget --no-check-certificate -c -O ${TEMPLATE_INPUT_FP} ${TEMPLATE_INPUT_URI}
export TEMPLATE_INPUT_BAND='1' # `ncks` makes dim=TSTEP first
# This script generates some basic statistics for a netCDF datavar.
STAT_SCRIPT_URI='https://github.com/TomRoche/GEIA_to_netCDF/raw/master/netCDF.stats.to.stdout.r'
STAT_SCRIPT_FN="$(basename ${STAT_SCRIPT_URI})"
export STAT_SCRIPT_FP="${TEST_DIR}/${STAT_SCRIPT_FN}"
wget --no-check-certificate -c -O ${STAT_SCRIPT_FP} ${STAT_SCRIPT_URI}
ls -alh ${TEST_DIR}/
R
# following done in R console-----------------------------------------
library(ncdf4)
library(raster)
library(maptools)
# run if rgeos not available
# gpclibPermit()
library(rgdal)
data(wrld_simpl) # from maptools
pdf.er <- Sys.getenv('PDF_VIEWER')
test.dir <- Sys.getenv('TEST_DIR')
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.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x=_0=-2556 +y_0=-1728'
out.fp <- Sys.getenv('DATA_OUTPUT_FP')
stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
source(stat.script.fp)
system(sprintf('ncdump -h %s', in.fp))
# netcdf GEIA_N2O_oceanic {
# 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) ;
# emi_n2o:units = "ton N2O-N/yr" ;
# emi_n2o:missing_value = -999. ;
# emi_n2o:long_name = "N2O emissions" ;
# emi_n2o:_FillValue = -999.f ;
# Note missing/fill values: more below.
netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=data.var.name)
# For /tmp/projectRasterTest/GEIA_N2O_oceanic.nc var=emi_n2o
# cells=64800
# obs=36143
# min=5.96e-08
# q1=30.4
# med=67.7
# mean=99.5
# q3=140
# max=1.17e+03
# Note min and max of the input, and that min > 0.
# make input raster
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
# make template raster
# truncate 'ncdump' to avoid IOAPI cruft
system(sprintf('ncdump -h %s | head -n 13', template.in.fp))
# netcdf emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT {
# dimensions:
# TSTEP = UNLIMITED ; // (25 currently)
# LAY = 1 ;
# ROW = 299 ;
# COL = 459 ;
# variables:
# float emi_n2o(TSTEP, LAY, ROW, COL) ;
# emi_n2o:long_name = "XYL " ;
# emi_n2o:units = "moles/s " ;
# emi_n2o:var_desc = "Model species XYL " ;
template.in.raster <- raster(template.in.fp, varname=data.var.name, band=template.band)
template.raster <- projectExtent(template.in.raster, out.crs)
# Warning message:
# In projectExtent(template.in.raster, out.crs) :
# 158 projected point(s) not finite
# is that "projected point(s) not finite" warning important? Probably not.
template.raster
# class : RasterLayer
# dimensions : 299, 459, 137241 (nrow, ncol, ncell)
# resolution : 53369.55, 56883.69 (x, y)
# extent : -12246449, 12250173, -4532510, 12475715 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +y_0=-1728 +ellps=WGS84
# values : none
# That looks good! Much more like the domain specification @
# https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
out.raster <-
projectRaster(
from=in.raster, to=template.raster, method='bilinear',
overwrite=TRUE, progress='window', format='CDF',
# args from ?writeRaster
# NAvalue=-999.0, # match emi_n2o:missing_value,_FillValue (TODO: copy)
NAflag=-999.0, # match emi_n2o:missing_value,_FillValue (TODO: copy)
varname=data.var.name,
varunit='ton N2O-N/yr',
longname='N2O emissions',
xname='COL',
yname='ROW',
filename=out.fp)
system(sprintf('ls -alht %s', test.dir))
system(sprintf('ncdump -h %s', out.fp))
# netcdf GEIA_N2O_oceanic_regrid {
# dimensions:
# COL = 459 ;
# ROW = 299 ;
# variables:
# double COL(COL) ;
# COL:units = "meter" ;
# double ROW(ROW) ;
# ROW:units = "meter" ;
# float emi_n2o(ROW, COL) ;
# emi_n2o:units = "ton N2O-N/yr" ;
# emi_n2o:missing_value = -3.4e+38f ;
# emi_n2o:_FillValue = -3.4e+38f ;
# emi_n2o:long_name = "N2O emissions" ;
# emi_n2o:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +y_0=-1728 +ellps=WGS84" ;
# emi_n2o:projection_format = "PROJ.4" ;
# emi_n2o:min = 5.9605e-08f ;
# emi_n2o:max = 783.92f ;
netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=data.var.name)
# For /tmp/projectRasterTest/GEIA_N2O_oceanic_regrid.nc var=emi_n2o
# cells=137241
# obs=53873
# min=5.96e-08
# q1=25.1
# med=55.5
# mean=81.1
# q3=93.1
# max=784
# plot to PDF
# do these once
date.format <- '%Y%m%d_%H%M'
pdf.dir <- test.dir
system(sprintf('mkdir -p %s', pdf.dir))
# repeat for each plot file
pdf.fn <- sprintf('test_%s.pdf', system(sprintf('date +%s', date.format), intern=TRUE))
pdf.fp <- sprintf('%s/%s', pdf.dir, pdf.fn)
pdf(file=pdf.fp, width=5.5, height=4.25)
# plot(out.raster, data.var.name) # only if RasterBrick or RasterStack
plot(out.raster)
# add a CONUS map
x <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]
us <- spTransform(x, CRS(out.crs))
plot(us, add=TRUE)
# close the device before using plotfile!
dev.off()
system(sprintf('ls -alh %s', pdf.fp))
# show the plot PDF
system(sprintf('%s %s &', pdf.er, pdf.fp))
# output is unexpectedly distorted
# try setting plot extents
pdf.fn <- sprintf('test_%s.pdf', system(sprintf('date +%s', date.format), intern=TRUE))
pdf.fp <- sprintf('%s/%s', pdf.dir, pdf.fn)
pdf(file=pdf.fp, width=5.5, height=4.25)
plot(out.raster, ext=template.raster)
# add a CONUS map
x <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]
us <- spTransform(x, CRS(out.crs))
plot(us, add=TRUE)
dev.off()
system(sprintf('ls -alh %s', pdf.fp))
system(sprintf('%s %s &', pdf.er, pdf.fp))
# no change from previous
# end example---------------------------------------------------------
The output I get is
http://tinyurl.com/rasterPlotExample_20121016
which
+ shows the US map, as expected
+ shows data around, but not inside, the map extents. Since this is
marine data, that's also good.
- shows a lot more than the expected map extents, i.e., a fringe
around CONUS
- shows the CONUS map very compressed, in a very small portion of the
plot extents
- has very unexpected legends, since the raster is col=459 x row=299.
How to get what I seek, i.e., with extents and shape more like
https://github.com/TomRoche/cornbeltN2O/wiki/images/CMAQ_AQMEII_MSFs.png
?
your assistance is appreciated, Tom Roche <Tom_Roche at pobox.com>
More information about the R-sig-Geo
mailing list