WCS from python

Basic interaction with the Web Coverage Service

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

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

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

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

print([op.name for op in wcs.operations])

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

print(list(wcs.contents))

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) 

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)

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

bdod_5_15_median.supportedFormats

bdod_5_15_median.boundingboxes

Obtaining a map segment from a WCS WCS version 1.0.0

Obtaining a map segment is another operation greatly simplified by the OWSLib package. This case is conducted from a pre-established bounding box.

A bounding box broadly matching Senegal:

bbox = (-1784000, 1356000, -1140000, 1863000)

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='bdod_5-15cm_median', 
    crs='urn:ogc:def:crs:EPSG::152160',
    bbox=bbox, 
    resx=250, resy=250, 
    format='GEOTIFF_INT16')

Now fetch the coverage for Senegal and save it to disk:

with open('./Senegal_bdod_5-15_median.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
ph = rasterio.open("Senegal_bdod_5-15_median.tif") 

Finally, use the plot class to plot the map:

from rasterio import plot
plot.show(ph, title='Mean pH between 0 and 5 cm deep in Senegal', cmap='gist_ncar')

Working with WCS 2.0

The WCS standard underwent a major update in 2017 that introduced important improvements. Among these is a more abstract treatment of data axes, thinking especially of space-time datasets and data cubes. The friendly OWSLib package has been updated accordingly and can be easily be used with the version 2.0 of the standard.

Connection set up and service inspection functions as before:

from owslib.wcs import WebCoverageService

wcs = WebCoverageService('http://maps.isric.org/mapserv?map=/map/phh2o.map',
                         version='2.0.1')

With the 2.0 version, supported formats are announced in slightly different way, let us thus inspect these in first place:

cov_id = 'phh2o_0-5cm_mean'
ph_0_5 = wcs.contents[cov_id]
ph_0_5.supportedFormats 

Perhaps the most noteciable change is the way segmentation is expressed. The infamous bounding boxes have been dropped at last and replaced by axis specific segmentation. The new argument subsets takes in a list of tuples, one tuple per axis. Each tuple is a triplet: (axis identifier, minimum, maximum). Here is an example with the same extent for Senegal used in the previous notebook:

subsets = [('X', -1784000, -1140000), ('Y', 1356000, 1863000)]

CRSs are also expressed in a different way in version 2.0, now referring to the opengis.net registry. Again the surrogate EPSG code 152160 is used, as the Petroleum folk remain unimpressed by the Homolosine projection

crs = "http://www.opengis.net/def/crs/EPSG/0/152160"

All information is now in place to issue a GetCoverage request and save the result to disk:

response = wcs.getCoverage(
  identifier=[cov_id], 
  crs=crs,
  subsets=subsets, 
  resx=250, resy=250, 
  format=ph_0_5.supportedFormats[0])

with open('./Senegal_pH_0-5_mean.tif', 'wb') as file:
  file.write(response.read())

The rasterio package can be used exactly as before to plot the segment obtained from the service:

import rasterio
from rasterio import plot

ph = rasterio.open("./Senegal_pH_0-5_mean.tif", driver="GTiff")
plot.show(ph, title='Mean pH between 0 and 5 cm deep in Senegal', cmap='gist_ncar')    

Last updated on 2025-10-16