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/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:
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:
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:
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:
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:
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:
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:
bdod_5_15_median = wcs.contents['bdod_5-15cm_Q0.5']
bdod_5_15_median.supportedCRS
[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]
bdod_5_15_median.supportedFormats
['GEOTIFF_INT16']
bdod_5_15_median.boundingboxes
[{'nativeSrs': 'EPSG:4326', 'bbox': (-179.991347553068, -55.9773009202418, 179.994461880094, 82.7192840534453)}, {'nativeSrs': 'EPSG:152160', 'bbox': (-19949000.0, -6147500.0, 19861750.0, 8361000.0)}]
Retrieving parameter for an area of interest
bbox = (32.25, -21.0, 35.0, -18.5) #coordinates given in latitude and longitude
wcs = WebCoverageService('http://maps.isric.org/mapserv?map=/map/sand.map', version='1.0.0')
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.
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:
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:
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:
from rasterio import plot
plot.show(sand_0_5, title='Mean Sand fraction between 0 and 5 cm deep', cmap='gist_ncar')
<Axes: title={'center': 'Mean Sand fraction between 0 and 5 cm deep'}>
from osgeo import gdal, osr
import matplotlib.pyplot as plt
import numpy as np
dataset = gdal.Open('./data/aoi_sand_0-5_mean.tif')
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
data_array = dataset.ReadAsArray().astype(np.int16)
print(data_array.dtype)
int16
Repeat this code to retrieve additional soil properties
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 (%)
Repeat this for clay, silt and bulk density for the top 5 cm soil layer