geoknife
has a number of output formats, but a new one
that was added in version 1.1.0 of the package (after the initial
release to CRAN) is the ability to output results as a zip file that
contains a series of geotiffs for each timestep requested. This vignette
is a brief introduction on how to do this using a few additional
packages that are in the Suggests
field in the
description.
see the geoknife
getting started vignette for more
details, but for the purposes of this vignette, we need to get some data
first and then plot it up.
First things first, load up the geoknife
package
library(geoknife)
For this example, I am going to use a dataset that is hosted by the
USGS, and includes an estimate of actual evapotranspiration. There is a
lot of data here, but I am just going to pluck a subset in time (one
timepoint for the month of July of one year) and a fairly large spatial
subset. I am going to use the OUTPUT_TYPE='geotiff'
argument that was added in version 1.1.0 of geoknife
.
<- webprocess(algorithm = list('OPeNDAP Subset' = "gov.usgs.cida.gdp.wps.algorithm.FeatureCoverageOPeNDAPIntersectionAlgorithm"))
knife
<- webdata(url='https://cida.usgs.gov/thredds/dodsC/ssebopeta/monthly',
fabric variable="et", #Monthly evapotranspiration
times=c('2014-07-15','2014-07-15'))
<- simplegeom(data.frame('point1' = c(-76,49), 'point2' = c(-93,40)))
stencil <- geoknife(stencil, fabric, knife, wait = TRUE, OUTPUT_TYPE = "geotiff") job
now that the job is complete (because wait=TRUE
was
used), we can download the result and unzip it:
<- file.path(tempdir(), 'et_data.zip')
dest
<- download(job, destination = dest, overwrite = TRUE)
file
unzip(file, exdir=file.path(tempdir(),'et'))
<- file.path(tempdir(),'et') tiff.dir
using the rasterVis
package, read this in and create a
raster object:
library(rasterVis)
library(raster)
<- raster(file.path(tiff.dir , dir(tiff.dir)))
et NAvalue(et) <- -1
Now, use ggplot2
and ggmap
to plot this as
a map:
library(ggmap)
library(ggplot2)
gplot(et, maxpixels = 5e5) +
geom_tile(aes(fill = value), alpha=1) +
scale_fill_gradientn("Actual ET (mm)", colours=rev(heat.colors(5))) +
coord_equal() +
theme_classic() +
theme(axis.line = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.ticks = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank()) +
scale_y_continuous(limits=c(et@extent@ymin, et@extent@ymax)) +
scale_x_continuous(limits=c(et@extent@xmin, et@extent@xmax))