# Access WCS with R

This tutorial will show how to access SoilGrids using Web Coverage Services and the R software. This tutorial will use the `XML` `terra` and `gdalUtilities` packages.

## Load libraries

```{r}
#| warning: false
#| eval: false

library(XML)
library(gdalUtilities)
library(terra)
```

## Set variables of interest

We define the variables for the soil property and layer of interest. See [here](SoilGrids_faqs.qmd#sec-faq2) for the naming conventions.

```{r}
#| warning: false
#| eval: false

voi = "nitrogen" # variable of interest
depth = "5-15cm"
quantile = "Q0.5"

voi_layer = paste(voi,depth,quantile, sep="_") # layer of interest 
```

We set other variables necessary for the WCS call for all kinds of requests.

```{r}
#| warning: false
#| eval: false

wcs_path = paste0("https://maps.isric.org/mapserv?map=/map/",voi,".map") # Path to the WCS. See maps.isric.org
wcs_service = "SERVICE=WCS"
wcs_version = "VERSION=2.0.1" # This works for gdal >=2.3; "VERSION=1.1.1" works with gdal < 2.3.
```

**NOTE**\
`wcs_version = "VERSION=2.0.1"` works for gdal \>=2.3; `wcs_version = "VERSION=1.1.1"` works with gdal \< 2.3.

## Example 1: Describe the coverage layer

###### 1. First we define the request as `DescribeCoverage` and we create a string for the full request using also the variables previously defined

```{r}
#| warning: false
#| eval: false

wcs_request = "request=DescribeCoverage" 

wcs = paste(wcs_path, wcs_service, wcs_version, wcs_request, sep="&")
```

###### 2. Then we create a XML that can be used with the `gdalinfo` utility after being saved to disk:

```{r}
#| warning: false
#| eval: false

l1 <- newXMLNode("WCS_GDAL")
l1.s <- newXMLNode("ServiceURL", wcs, parent=l1)
l1.l <- newXMLNode("CoverageName", voi_layer, parent=l1)

# Save to local disk
xml.out = "./sg.xml"
saveXML(l1, file = xml.out)
```

###### 3. Finally we use `gdalinfo` to get the description of the layer

```{r}
#| warning: false
#| eval: false

gdalinfo("./sg.xml")
```

## Example 2: Download a Tiff for a small region of interest (ROI)

The terra package can be used to directly download the Tiff for a small ROI, however for larger regions it throws the error `Error: [rast] file does not exist:..`. If this is the case take a look at Example 3 for how to download a larger Tiff for a larger ROI.

A couple of additional arguments need to be added to the url and the request type needs to be changed from `DescribeCoverage`to `GetCoverage`. The format (GEOTIFF_INT16) and the region of interest defined need to be defined.

```{r}
#| warning: false
#| eval: false

wcs_request = "request=GetCoverage"
format = "format=GEOTIFF_INT16"
roi = "subset=X(-1026935.24156938,-1015591.34251064)&subset=Y(2005364.48088881,2014307.50831563)"
layer = paste0("coverageid=",voi_layer)
url_request = paste(wcs_path, layer, wcs_service, wcs_version, wcs_request,format,roi, sep="&" )
```

The url we have just defined can then be read using the rast function of the terra package. Our Tiff can then be viewed and saved as needed.

```{r}
#| warning: false
#| eval: false

small_tiff = rast(url_request)

plot(small_tiff)

writeRaster(small_tiff, './rast.tif')
```

## Example 3: Download a Tiff for a small region of interest (ROI)

###### 1. First we define the request with the parameters for the region of interest, projection, format, resolution

An example (Ghana) using Homolosine projection.

```{r}
#| warning: false
#| eval: false

bb=c(-337500.000,1242500.000,152500.000,527500.000) # Example bounding box (homolosine)
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' # proj string for Homolosine projection
```

We paste all the components of the WCS request:

```{r}
#| warning: false
#| eval: false

wcs = paste(wcs_path,wcs_service,wcs_version,sep="&") # This works for gdal >= 2.3
```

**NOTE** For gdal \< 2.3 an extra `&` is necessary at the end of the string: `wcs = paste0(wcs,"&")`

###### 2. Then we create a XML that can be used with the `gdal` utility after being saved to disk:

```{r}
#| warning: false
#| eval: false

l1 <- newXMLNode("WCS_GDAL")
l1.s <- newXMLNode("ServiceURL", wcs, parent=l1)
l1.l <- newXMLNode("CoverageName", "nitrogen_5-15cm_Q0.5", parent=l1)

# Save to local disk
xml.out = "./sg.xml"
saveXML(l1, file = xml.out)
```

###### 3. Finally we use `gdal_translate` to get the geotiff locally.

```{r}
#| warning: false
#| eval: false

# Download raster as GeoTIFF (Warning: it can be large!)
file.out <- './test.tif'

gdal_translate(xml.out, file.out,
    tr=c(250,250), projwin=bb,
    projwin_srs =igh, co=c("TILED=YES","COMPRESS=DEFLATE","PREDICTOR=2","BIGTIFF=YES")
)
```

The downloaded file can then be read into R.

See [this page](WebDav.qmd) for examples on how to reproject.

# Common errors:

If you get an SSL error (for example: `SSL certificate problem: self signed certificate in certificate chain`) please try to add the following config parameter to the gdal call: `--config GDAL_HTTP_UNSAFESSL YES`:

```{r}
#| warning: false
#| eval: false

gdal_translate(xml.out, file.out,
    tr=c(250,250), projwin=bb,
    projwin_srs =igh, co=c("TILED=YES","COMPRESS=DEFLATE","PREDICTOR=2","BIGTIFF=YES","GDAL_HTTP_UNSAFESSL=YES")
)
```
