# WCS from python 

## Basic interaction with the Web Coverage Service

The interaction with a WCS is made very convenient by the [OWSLib](https://owslib.readthedocs.io/en/latest/index.html) 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:


```python
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:

```python
print([op.name for op in wcs.operations])
```
The full list of coverages available from this service is in the `contents` property:

```python
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:

```python
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:

```python
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](https://www.globalsoilmap.net/) specifications.

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

```python
bdod_5_15_median = wcs.contents['bdod_5-15cm_Q0.5']
bdod_5_15_median.supportedCRS
```

```python
bdod_5_15_median.supportedFormats
```

```python
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:

```python
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.

```python
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:

```python
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:

```python
import rasterio
ph = rasterio.open("Senegal_bdod_5-15_median.tif") 
```

Finally, use the plot class to plot the map:

```python
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](https://geopython.github.io/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:

```python
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:
 
```python
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:

```python
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](http://www.opengis.net/) registry. Again the surrogate EPSG code 152160 is used, as the Petroleum folk remain unimpressed by the Homolosine projection
 
```python
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:
 
```python
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:
 
```python    
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')    
```

