Basic interaction with the Web Coverage Service¶

https://soilgrids.readthedocs.io/en/latest/

https://www.isric.org/explore/soilgrids/faq-soilgrids#What_is_SoilGrids

https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/01-WCS-basics.ipynb?ref_type=heads

https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/markdown/access_on_gee.md

The interaction with a WCS is made very convinient by the OWSLib pacakge. This section shows how to use OWSLib to obtain relevant information about the service and the maps it serves.

Available prediction layers:

Capture.JPG

Part 1¶

First load the WebCoverageService class from the OWSLib and create a connection to a service, in this case the one serving the predictions for sand fraction:

In [1]:
from owslib.wcs import WebCoverageService
wcs = WebCoverageService('http://maps.isric.org/mapserv?map=/map/bdod.map', version='1.0.0')

The absense of errors means that a connection was sucessfully established. To get more information about this service, start by identifying the operations available:

In [2]:
print([op.name for op in wcs.operations])
['GetCapabilities', 'DescribeCoverage', 'GetCoverage']

The full list of coverages available from this service is in the contents property:

In [3]:
print(list(wcs.contents))
['bdod_0-5cm_Q0.5', 'bdod_0-5cm_Q0.05', 'bdod_0-5cm_Q0.95', 'bdod_0-5cm_mean', 'bdod_0-5cm_uncertainty', 'bdod_5-15cm_Q0.5', 'bdod_5-15cm_Q0.05', 'bdod_5-15cm_Q0.95', 'bdod_5-15cm_mean', 'bdod_5-15cm_uncertainty', 'bdod_15-30cm_Q0.5', 'bdod_15-30cm_Q0.05', 'bdod_15-30cm_Q0.95', 'bdod_15-30cm_mean', 'bdod_15-30cm_uncertainty', 'bdod_30-60cm_Q0.05', 'bdod_30-60cm_Q0.5', 'bdod_30-60cm_Q0.95', 'bdod_30-60cm_mean', 'bdod_30-60cm_uncertainty', 'bdod_60-100cm_Q0.05', 'bdod_60-100cm_Q0.5', 'bdod_60-100cm_Q0.95', 'bdod_60-100cm_mean', 'bdod_60-100cm_uncertainty', 'bdod_100-200cm_Q0.05', 'bdod_100-200cm_Q0.5', 'bdod_100-200cm_Q0.95', 'bdod_100-200cm_mean', 'bdod_100-200cm_uncertainty']

That is a large set of coverages, but it is easy to filter the dictionary. For instance, to get the name of all coverages for the 0 cm to 5 cm depth interval:

In [4]:
names = [k for k in wcs.contents.keys() if k.startswith("bdod_0-5cm")]
print(names)
['bdod_0-5cm_Q0.5', 'bdod_0-5cm_Q0.05', 'bdod_0-5cm_Q0.95', 'bdod_0-5cm_mean', 'bdod_0-5cm_uncertainty']

Or to search for all the coverages reporting the median prediction:

In [5]:
q0_5_covs = [k for k in wcs.contents.keys() if k.find("Q0.5") != -1]
print(q0_5_covs)
['bdod_0-5cm_Q0.5', 'bdod_5-15cm_Q0.5', 'bdod_15-30cm_Q0.5', 'bdod_30-60cm_Q0.5', 'bdod_60-100cm_Q0.5', 'bdod_100-200cm_Q0.5']

These are the SoilGrids predictions for bulk density for the six standard depths defined in the GlobalSoilMap specifications.

The details for one of these coverages can be inspected using the identifiers above:

In [6]:
bdod_5_15_median = wcs.contents['bdod_5-15cm_Q0.5']
bdod_5_15_median.supportedCRS
Out[6]:
[urn:ogc:def:crs:EPSG::152160,
 urn:ogc:def:crs:EPSG::4326,
 urn:ogc:def:crs:EPSG::3857,
 urn:ogc:def:crs:EPSG::54009,
 urn:ogc:def:crs:EPSG::54012,
 urn:ogc:def:crs:EPSG::152160]
In [7]:
bdod_5_15_median.supportedFormats
Out[7]:
['GEOTIFF_INT16']
In [8]:
bdod_5_15_median.boundingboxes
Out[8]:
[{'nativeSrs': 'EPSG:4326',
  'bbox': (-179.991347553068,
   -55.9773009202418,
   179.994461880094,
   82.7192840534453)},
 {'nativeSrs': 'EPSG:152160',
  'bbox': (-19949000.0, -6147500.0, 19861750.0, 8361000.0)}]

Part 2¶

Retrieving parameter for an area of interest

In [9]:
bbox = (32.25, -21.0, 35.0, -18.5) #coordinates given in latitude and longitude
In [10]:
wcs = WebCoverageService('http://maps.isric.org/mapserv?map=/map/sand.map', version='1.0.0')
In [11]:
q0_5_covs = [k for k in wcs.contents.keys() if k.find("Q0.5") != -1]
print(q0_5_covs)
['sand_0-5cm_Q0.5', 'sand_5-15cm_Q0.5', 'sand_15-30cm_Q0.5', 'sand_30-60cm_Q0.5', 'sand_60-100cm_Q0.5', 'sand_100-200cm_Q0.5']

The getCoverage method can now be used to fetch the map segment within the bounding box. Note the other parameters, Section 1 showed how to obtain them.

In [12]:
response = wcs.getCoverage(
    identifier='sand_0-5cm_Q0.5',
    crs='urn:ogc:def:crs:EPSG::4326', 
    bbox=bbox, 
    resx=0.002, resy=0.002, 
    format='GEOTIFF_INT16')

Now fetch the coverage for Area of Interest (as specified by bbox) and save it to disk:

In [13]:
with open('./data/aoi_sand_0-5_mean.tif', 'wb') as file:
    file.write(response.read())

With the data on the client side some regular interaction can be started with a library like rasterIO. First open the file from disk:

In [14]:
import rasterio
sand_0_5 = rasterio.open("./data/aoi_sand_0-5_mean.tif", driver="GTiff")

Finally, use the plot class to plot the map:

In [15]:
from rasterio import plot
plot.show(sand_0_5, title='Mean Sand fraction between 0 and 5 cm deep', cmap='gist_ncar')
Out[15]:
<Axes: title={'center': 'Mean Sand fraction between 0 and 5 cm deep'}>

Now using GDAL, Numpy and Matplotlib¶

In [16]:
from osgeo import gdal, osr
import matplotlib.pyplot as plt
import numpy as np
In [17]:
dataset = gdal.Open('./data/aoi_sand_0-5_mean.tif')
In [18]:
coverage_info = gdal.Info(dataset)
print (coverage_info)
Driver: GTiff/GeoTIFF
Files: ./data/aoi_sand_0-5_mean.tif
Size is 1375, 1250
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (32.250000000000000,-18.500000000000000)
Pixel Size = (0.002000000000000,-0.002000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  32.2500000, -18.5000000) ( 32d15' 0.00"E, 18d30' 0.00"S)
Lower Left  (  32.2500000, -21.0000000) ( 32d15' 0.00"E, 21d 0' 0.00"S)
Upper Right (  35.0000000, -18.5000000) ( 35d 0' 0.00"E, 18d30' 0.00"S)
Lower Right (  35.0000000, -21.0000000) ( 35d 0' 0.00"E, 21d 0' 0.00"S)
Center      (  33.6250000, -19.7500000) ( 33d37'30.00"E, 19d45' 0.00"S)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray

In [19]:
data_array = dataset.ReadAsArray().astype(np.int16)
print(data_array.dtype)
int16

Repeat this code to retrieve additional soil properties

In [20]:
plt.figure(figsize=(8,8))
plt.imshow(data_array, interpolation = 'nearest', cmap = 'gist_ncar')#, vmin=0, vmax=20)
plt.title('Mean Sand fraction between 0 and 5 cm deep')
plt.colorbar(shrink=0.4)
#plt.grid(False)
plt.show()

# Note unit of Proportion of sand particles (> 0.05 mm) in the fine earth fraction
# Mapped units (pixel value) = g/kg, Conversion factor = 10, Conventional units = g/100g (%)

Exercise:¶

Repeat this for clay, silt and bulk density for the top 5 cm soil layer

In [ ]: