Downloading SoilGrids data globally using R

I’m trying to download the ISRIC SoilGrids data via R (Version 3.6.2 on macOS Catalina Version 10.15.5) as described in the tutorial here:
https://www.isric.org/explore/soilgrids/soilgrids-access

Following this guide I can download the data for the small example bounding box. However, I’m interested in the global data so my modified R script looks like this:

library(rgdal)
library(gdalUtils)

bb=c(-20037500,-6729000,20037500,8600750) #bounding box parameters (in Homolosine) for whole globe, as suggested on the website 
igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' #proj string for Homolosine projection


gdal_translate(
  '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt',
  "./crop_roi_igh_r.vrt",
  of="VRT",
  tr=c(25000,25000),
  projwin=bb,
  projwin_srs =igh,
  verbose=TRUE
)


gdalwarp(
  "./crop_roi_igh_r.vrt",
  "./crop_roi_ll_r.vrt",
  s_src=igh, 
  t_srs="EPSG:4326", 
  of="VRT"
)

gdal_translate(
  "./crop_roi_ll_r.vrt",  
  "./crop_roi_ll_r.tif", 
  co=c("TILED=YES","COMPRESS=DEFLATE","PREDICTOR=2","BIGTIFF=YES")
)

Running

gdal_translate(
  '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt',
  "./crop_roi_igh_r.vrt",
  of="VRT",
  tr=c(25000,25000),
  projwin=bb,
  projwin_srs =igh,
  verbose=TRUE
)

aborts with the following error message:

Checking gdal_installation…
GDAL version 2.4.2
GDAL command being used: "/usr/local/Cellar/gdal/2.4.2_3/bin/gdal_translate" -tr 2500 2500 -projwin -20037500 -6729000 20037500 8600750 -of "GTiff" -projwin_srs "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs" "/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt" "./crop_roi_igh_r.tif"
ERROR 1: -projwin_srs ignored since coordinate transformation failed.
Input file size is 159246, 58034

What I tried so far:

1. incrementally increasing the bounding box from bb=c(-337500.000,1242500.000,152500.000,527500.000)
   up to bb=c(-19861874.000, 1242500.000, 19861874.000, -1500000.000)

This spans a belt from -179.xx to 179.xx Longitude but I can’t increase Latitude parameters even to 20°N/S before encountering the error.

2. Going for a different parameter (mean clay content, which I’m actually more interested in) via

gdal_translate(
  '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/clay/clay_0-5cm_mean.vrt',
  "./crop_clay_igh_r.vrt",
  of="VRT",
  tr=c(25000,25000),
  projwin=bb,
  projwin_srs =igh,
  verbose=TRUE
)

This gives the error

ERROR 1: Error: Computed -srcwin 71784 36135 16030 -6132 has negative width and/or height.
Input file size is 159246, 58034