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

Merge branch 'wcsfromr' into 'master'

changed tutorial markdown/webdav_from_R.md

See merge request !1
parents fabb1154 1e868797
Branches
No related tags found
1 merge request!1changed tutorial markdown/webdav_from_R.md
...@@ -4,14 +4,13 @@ ...@@ -4,14 +4,13 @@
# [WebDAV](https://files.isric.org/soilgrids/data/recent/): direct access to the maps in VRT format. # [WebDAV](https://files.isric.org/soilgrids/data/recent/): direct access to the maps in VRT format.
### R ### R
**This code is tested for GDAL version 2.4.x. When running with GDAL version 3.x.x there may be problems. We are working on providing updates**
#### Load libraries, select boundary box and cell size #### Load libraries, select boundary box and cell size
Initially you need to load the following libraries and select the bounding box of the area you are interested in: Initially you need to load the following libraries and select the bounding box of the area you are interested in:
```R ```{r}
library(rgdal) library(terra)
library(gdalUtils) library(gdalUtilities)
bb=c(-337500.000,1242500.000,152500.000,527500.000) # Example bounding box (homolosine) for Ghana bb=c(-337500.000,1242500.000,152500.000,527500.000) # Example bounding box (homolosine) for Ghana
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' # proj string for Homolosine projection igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' # proj string for Homolosine projection
...@@ -21,13 +20,12 @@ sg_url="/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.o ...@@ -21,13 +20,12 @@ sg_url="/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.o
#### To a geotiff in Homolosine #### To a geotiff in Homolosine
This GDAL command will create a local geotiff in the Homolosine projection This GDAL command will create a local geotiff in the Homolosine projection
```R ```{r}
gdal_translate(paste0(sg_url,'ocs/ocs_0-30cm_mean.vrt'), gdal_translate(paste0(sg_url,'ocs/ocs_0-30cm_mean.vrt'),
"./crop_roi_igh_r.tif", "./crop_roi_igh_r.tif",
tr=c(250,250), tr=c(250,250),
projwin=bb, projwin=bb,
projwin_srs =igh, projwin_srs =igh)
verbose=TRUE)
``` ```
#### To a geotiff in a different projection #### To a geotiff in a different projection
...@@ -36,30 +34,30 @@ The following commands describe a workflow to obtain a VRT or a GeoTiff for an a ...@@ -36,30 +34,30 @@ The following commands describe a workflow to obtain a VRT or a GeoTiff for an a
##### To local VRT in homolosine (directly from the webdav connection) ##### To local VRT in homolosine (directly from the webdav connection)
The first step is to obtain a VRT for the area of interest in the Homolosine projection. We suggest to use VRT for the intermediate steps to save space and computation times. The first step is to obtain a VRT for the area of interest in the Homolosine projection. We suggest to use VRT for the intermediate steps to save space and computation times.
```R ```{r}
gdal_translate(paste0(sg_url,'ocs/ocs_0-30cm_mean.vrt'), gdal_translate(paste0(sg_url,'ocs/ocs_0-30cm_mean.vrt'),
"./crop_roi_igh_r.vrt", "./crop_roi_igh_r.vrt",
of="VRT",tr=c(250,250), of="VRT",tr=c(250,250),
projwin=bb, projwin=bb,
projwin_srs =igh, projwin_srs =igh)
verbose=TRUE)
``` ```
##### To a VRT in, for example, LatLong ##### To a VRT in, for example, LatLong
The following command will generate a VRT in the projection of your choice: The following command will generate a VRT in the projection of your choice:
```R ```{r}
gdalwarp("./crop_roi_igh_r.vrt", gdalwarp("./crop_roi_igh_r.vrt",
"./crop_roi_ll_r.vrt", "./crop_roi_ll_r.vrt",
s_src=igh, s_srs=igh,
t_srs="EPSG:4326", t_srs="EPSG:4326",
of="VRT") of="VRT",
overwrite = TRUE)
``` ```
##### To a final Geotiff ##### To a final Geotiff
The following command will generate a Geotiff in the projection of your choice for the area of interest defined above The following command will generate a Geotiff in the projection of your choice for the area of interest defined above
```R ```{r}
gdal_translate("./crop_roi_ll_r.vrt", gdal_translate("./crop_roi_ll_r.vrt",
"./crop_roi_ll_r.tif", "./crop_roi_ll_r.tif",
co=c("TILED=YES","COMPRESS=DEFLATE","PREDICTOR=2","BIGTIFF=YES")) co=c("TILED=YES","COMPRESS=DEFLATE","PREDICTOR=2","BIGTIFF=YES"))
...@@ -68,14 +66,18 @@ gdal_translate("./crop_roi_ll_r.vrt", ...@@ -68,14 +66,18 @@ gdal_translate("./crop_roi_ll_r.vrt",
#### Read in R #### Read in R
Finally you can read any of the generated VRTs or GeoTiffs in R for further analysis Finally you can read any of the generated VRTs or GeoTiffs in R for further analysis
```R ```{r}
r=readGDAL("./crop_roi_igh_r.vrt") # Or any other of the files produce above rst=rast("./crop_roi_igh_r.vrt") # Or any other of the files produce above
```
```{r}
plot(rst) # Plot the file
``` ```
# To download a global geotiff in Homolosine projection # To download a global geotiff in Homolosine projection
```R ```{r}
gdal_translate(paste0(sg_url,'ocs/ocs_0-30cm_mean.vrt'), gdal_translate(paste0(sg_url,'ocs/ocs_0-30cm_mean.vrt'),
"./crop_roi_igh_r.tif", "./crop_roi_igh_r.tif")
verbose=TRUE)
``` ```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment