# Quick example of accessing a WCS using gdal via R
# assumes gdal is available in PATH!
# Edward P. Morris (UCA), 2017-01-30.
# example adapted from: http://gsif.isric.org/doku.php/wiki:tutorial_soilgrids

# Load packages
require(XML)
require(rgdal)
require(raster)

# Define WCS XML file
wcs = "http://fast.openearth.eu/geoserver/ows?service=WCS&version=2.0.1&request=GetCapabilities"
l1 <- newXMLNode("WCS_GDAL")
l1.s <- newXMLNode("ServiceURL", wcs, parent=l1)
l1.l <- newXMLNode("CoverageName", "England_Study_Site_DonnaNook:DTM", parent=l1)
l1

# Save to local disk
xml.out = "~/Desktop/England_Study_Site_DonnaNook:DTM.xml"
saveXML(l1, file = xml.out)

# Check the XML file and WCS definition is correct
system(paste('gdalinfo', xml.out))

# Download raster as GeoTIFF (92 MB!)
file.out <- '~/Desktop/England_Study_Site_DonnaNook:DTM.tif'
system(paste('gdal_translate', ' ', xml.out, file.out))

# Plot
r <- raster('~/Desktop/England_Study_Site_DonnaNook:DTM.tif')
NAvalue(r) <- -9999
plot(r, xlab="Longitude", ylab="Latitude", main="Elevation (m)")