# WebDAV: R - terra

[WebDAV](https://files.isric.org/soilgrids/latest/data/): direct access to the maps using terra.

### Access SoilGrids 2.0 layers using terra

For access to the *World Refernce Base* (WRB) layers please skip to @sec-WRB.

By following the outlined steps, you will be able to access and download your desired SoilGrids layer using R. In the initial two steps, you will define the necessary variables. The subsequent steps (steps 3 and 4) utilise these defined variables and can be executed by the user without any need for modifications. In step 2 you will define a region of interest (roi). If you would like global data then skip to step 3b. *Please note that downloading global data may take quite some time. SoilGrids is intended for global applications, however in some cases it is appropriate to use a roi, i.e., continental scale analyses.*

#### Step 1 - Load library & select layer of interest

To begin, make sure you load the "terra" library, define the SoilGrids' layer of interest and desired resolution and specify where you wish to save your file. You can find a list of available layers at [**https://files.isric.org/soilgrids/latest/data/**](https://files.isric.org/soilgrids/latest/data/).

SoilGrids native resolution is 250m but aggregated products for the **mean only** of selected soil properties are available at 1000m and 5000m as .tif files , check the availability here: [**https://files.isric.org/soilgrids/latest/data_aggregated/**](https://files.isric.org/soilgrids/latest/data_aggregated/)

##### Load the terra package and define output paths

``` {r}
#| warning: false
#| eval: false
library(Rcpp)
library(terra)

# EDIT: Enter your desired resolution and output options
resolution = "250m" # or "1000m"" or "5000m"
output_folder = "output/folder/location/" # enter here your filepath to your output folder
output_filename = "output.tif" # choose a name for your output file
```

##### Define layer of interest for soil properties

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

# EDIT - the voi, depth and quantile
voi = "nitrogen" # variable of interest
depth = "5-15cm" # 0-5cm, 5-15cm, 15-30cm, 30-60cm, 60-100cm or 100-200cm
quantile = "mean" # mean, Q0.05, Q0.5, Q0.95 or uncertainty

# RUN - Here we define the variables for the layer name and coordinate reference system. 
voi_layer = paste(voi,depth,quantile, sep="_") # layer of interest
rst_crs = "ESRI:54052" # ESRI code of the interupted Goode Homolosine projection used by SoilGrids
```


#### Step 2 - Define Extent

Next, you will define the extent of the region of interest (ROI). You have two options for this:

1.  **Step 2a:** If you are already familiar with the minimum and maximum x and y coordinates, along with the Coordinate Reference System (CRS) of your ROI, you can manually input these values.

2.  **Step 2b:** Alternatively, you can extract the extent from a raster or vector file that represents your ROI.

##### 2a Get bounding box from known extent for region of interest

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

in_crs="EPSG:4326" # CRS of the bounding box. use EPSG code or WKT 

xmin = 0 #bounding box xmin
xmax = 2 #bounding box xmax
ymin = 40 #bounding box ymin
ymax = 42 #bounding box ymax

bb=ext(c(xmin,xmax,ymin,ymax))
```

##### 2b Get bounding box from raster or vector file of region of interest

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

roi_path = "roi/file/path/roi_filename.gpkg"

roi = vect(roi_path) #use vect(filename) for vector data or rast(filename) for raster

bb = ext(roi)

in_crs = crs(roi)
```

#### Step 4 - Download global SoilGrids layer 

In this step the SoilGrids layer is loaded in terra. The 'window' function is used to select the region of interest. 

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

if (resolution == "250m") { 
rstFile = paste0("/vsicurl/https://files.isric.org/soilgrids/latest/data/", voi, '/', voi_layer,'.vrt')
} else if (resolution == "1000m") {
rstFile = paste0("/vsicurl/https://files.isric.org/soilgrids/latest/data_aggregated/1000m/", voi, '/', voi_layer,'_1000.tif')
} else {
rstFile = paste0("/vsicurl/https://files.isric.org/soilgrids/latest/data_aggregated/5000m/", voi, '/', voi_layer,'_5000.tif')
}

rst = rast(rstFile)
crs(rst) = rst_crs # assign crs as the ESRI code rather than proj string

plot(rst) # view layer of interest for the whole world (may take a minute or ten)
```

#### Step 4 - Subset SoilGrids layer of interest to given ROI

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

bb_proj = project(bb,from=crs(in_crs),to=crs(rst_crs)) #project bb to same crs as raster layer

window(rst) = bb_proj # get just roi

plot(rst) # quick visual check 
```

#### Step 4 - Save file

Finally, you can optionally re-project the layer to the CRS of the ROI (un-comment the line if required), if not, go ahead and save your file.

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

# project to initial crs of roi                    
#rst = project(rst, in_crs) #optional

# save  your file                  
writeRaster(rst, paste0(output_folder, output_filename), overwrite = TRUE)
```


## Access the WRB layers using terra {#sec-WRB}

##### 1b Define layer of interest for soil class

In SoilGrids 2.0 the Soil type classes were not re done XXX what is this sentance. (For more information see the relevant scetions of the FAQs here XXX)
The WRB layers differ to the SoilGrids 2.0 layers in the WEBDAV as they are only available at 250m and have the crs code of epsg:4236.
The process to acess them is similar to that as above but with some important tweaks. 

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

# EDIT - enter you soil class of choice or 'MostProbable'
soil_class = 'Acrisols' # choose soil class or 'MostProbable'

# RUN - We define the inputs (depth and quantile are not relevant here)
voi = 'wrb' # this should not change
voi_layer = soil_class 
rst_crs = 'EPSG:4326' # EPSG code for the WRB layers
```
