Skip to content
Snippets Groups Projects
Commit 23fe857b authored by Genova, Giulio's avatar Genova, Giulio
Browse files

adding markdown/webdav_from_R_terra.md

parent f091a51a
No related branches found
No related tags found
No related merge requests found
<a href="https://www.isric.org" rel="isric.org"> <img src="https://www.isric.org/themes/custom/basic/logo.svg" height="130"/> <a href="https://soilgrids.org" rel="soilgrids.org"> <img src="https://www.isric.org/sites/default/files/styles/gallery_big_image_900x700/public/SoilGrids_banner_web.png" height="130"/>
# [WebDAV](https://files.isric.org/soilgrids/data/recent/): direct access to the maps using terra.
### R
By following the four steps outlined below, you will be able to access and download your desired SoilGrids layer using the terra package in R. Step 1 involves specifying the output file path and defining necessary variables such as resolution, variable of interest, depth, and quantile to download your selected SoilGrids layer. In Step 2 you will define a region of interest (roi)\*. If you require global data then skip Step 2 and 3b. Steps 3 and 4 utilise the variables defined in Step 1 (& optionally Step 2) to download your SoilGrids layer. Steps 3 & 4 can be executed by the user without any need for modifications.
*\*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, load the "terra" library, define the required resolution, output location and output filename, and specify the SoilGrids' layer of interest. 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/)
If your variable of interest is a **soil property** then this should be defined in **Step 1a.** If your variable of interest is a soil class (reference soil group) from the **World Reference Base** (WRB) then this should be defined in **Step 1b**. Please note that SoilGrids and WRB use different coordinate-reference systems, take care that the correct one is assigned ("ESRI:54052" for soil properties; 'EPSG:4326' for WRB reference groups).
``` r
library(terra)
resolution = "250m" # or "1000m"" or "5000m"
output_folder = "output/folder/location/" # filepath to output folder on your computer
output_filename = "output.tif" # choose a name for your file
```
##### **1a Define layer of interest for soil properties**
``` r
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
voi_layer = paste(voi,depth,quantile, sep="_") # layer of interest
rst_crs = "ESRI:54052" # ESRI code of the interrupted Goode Homolosine projection used by SoilGrids
```
##### **1b Define layer of interest for soil class**
*Only available at 250m resolution, please make sure resolution is 250m in Step 1.*
``` r
voi = "wrb" # this should not change
soil_class = "Acrisols" # choose soil class or "MostProbable"
voi_layer = soil_class
rst_crs = "EPSG:4326" # EPSG code for the WRB layers
```
#### 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 (longitude) and y (latitude) 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
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
roi_path = "roi/file/path/roi_filename.gpkg" #set the file path to your roi file (vector or raster format)
roi = vect(roi_path) # use vect(filename) for vector data or rast(filename) for raster
bb = ext(roi)
in_crs = crs(roi)
```
#### Step 3 - Download SoilGrids layer
In this step the SoilGrids layer is loaded in terra. In step 3b the 'window' function is used to select the region of interest. Skip step 3b if a global coverage is required.
##### 3a Get global SoilGrids layer of interest
``` r
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)
```
##### 3b Get SoilGrids layer of interest for given ROI
``` r
bb_proj = project(bb,from=crs(in_crs),to=crs(rst_crs)) #project bb to same crs as raster layer
window(rst) = NULL # remove any existing window
window(rst) = bb_proj # get just roi
plot(rst) # quick visual check
```
#### Step 4 - Save file
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
# project to initial crs of roi
#rst = project(rst, in_crs) #optional
# save your file
writeRaster(rst, paste0(output_folder, output_filename), overwrite = TRUE)
```
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment