Reading VRT from SoilGrids in R

Question

I need the global SOC stocks map at 0-30cm from the new SoilGrids 2.0: https://files.isric.org/soilgrids/latest/data/ocs/

I have no previous experience with VRT files. I have simply downloaded "ocs_0-30cm_mean.vrt" and "ocs_0-30cm_mean.vrt.ovr". Honestly, I don't know if I need one, the other, or both.

In R, I load it with the raster package.

soc <- raster("ocs_0-30cm_mean.vrt.ovr") # t/ha
soc
class      : RasterLayer 
dimensions : 14509, 39812, 577632308  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 39812, 0, 14509  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : ocs_0-30cm_mean.vrt.ovr 
names      : ocs_0.30cm_mean.vrt 
values     : -32768, 32767  (min, max)

However, the resolution appears to be 1 degree, when it should be 250m. The extent makes no sense (I guess it should be c(-180, 180, -90, 90), and there is no CRS.

I can plot (soc), though I don't know how to change the resolution to e.g. 0.25. I have also tried the same with "ocs_0-30cm_mean.vrt" (instead of "ocs_0-30cm_mean.vrt.ovr"), but the file does not plot

plot(soc)
Error in rgdal::getRasterData(con, offset = offs, region.dim = c(1, nc),  : 
  Failure during raster IO

How do I make this raster work, with 0.25 resolution and appropriate extent?

Answer 1

It is better that you consult the Soilgrids FAQ mainly the section How can I download Soilgrids the webdav download describes what are the vrt and ovr files.

vrt is a virtual XML file that points to tiles creating a mosaic that behaves like a single file, therefore if you copy only the vrt to your computer you will not get anything, ovr are overviews of the vrt file

The strange numbers are the homolosine projection give it a look on the FAQ Homolosine

I am not an R programmer but you can access the files on server as a GDAL virtual file system (as if the file would be on your computer and in this case it is /vsicurl

link<-"/vsicurl/https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt"
> raster(link)
class      : RasterLayer 
dimensions : 58034, 159246, 9241682364  (nrow, ncol, ncell)
resolution : 250, 250  (x, y)
extent     : -19949750, 19861750, -6147500, 8361000  (xmin, xmax, ymin, ymax)
crs        : +proj=igh +ellps=WGS84 +units=m +no_defs 
source     : /vsicurl/https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt 
names      : ocs_0.30cm_mean 
values     : 6, 208  (min, max)

Pay attention that each layer is 5 gigas compressed

---

Using gdalwarp it is possible to convert the remote vrt file into a local geotiff.

gdalwarp -t_srs EPSG:4326 -multi -wm 200 -co BIGTIFF=YES -co COMPRESS=DEFLATE -co TILED=TRUE \
"/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt"
"soc_0-5cm_mean.tif"

Just take create that compression is active: -co COMPRESS=DEFLATE

Answer 2

Based on the answer from Jorge Mendes, I can use gdalUtils in R to read the online VRT file, transform it to a more friendly projection and resolution, and store a copy in my local disk as a .tif file:

# https://files.isric.org/soilgrids/latest/data/ocs/
library(gdalUtils)
gdalwarp(t_srs="EPSG:4326", multi=TRUE, wm=200, 
co=c("BIGTIFF=YES", "COMPRESS=DEFLATE", "TILED=TRUE"),
tr=c(0.25,0.25), # Desired output resolution
verbose=T,
"/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/soc/soc_0-5cm_mean.vrt", # Input VRT
"soc_0-5cm_SoilGrids2.tif") # Output file