We will use the python libraries listed below. They can be installed with following command pip install os owslib gdal requests pandas matplotlib basemap.

import os
from owslib.wcs import WebCoverageService
from osgeo import gdal
import requests
import xml.dom.minidom
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

1 Working with data from EMODnet bathymetry

1.1 Access raster data through WCS

We use the owslib package to access the EMODnet bathymetry WCS.

# define the connection
url = "https://ows.emodnet-bathymetry.eu/wcs?"
wcs = WebCoverageService(url, version='1.0.0', timeout = 320)
print(wcs.identification.type)
## OGC:WCS
print(wcs.identification.title)
## EMODnet Bathymetry WCS

We set the desired parameters of the bathymetry raster we want to download.

# define variables
clipfile =  r'temp.tif'
requestbbox = (2,51.5,5,54)
layer = 'emodnet:mean'
# get the data
bathyml = 'emodnet:mean'
sed = wcs[layer] # this is necesaary to get essential metadata from the advertised layer
print(sed.keywords)
## ['emodnet_v9_mean', 'WCS', 'GeoTIFF']
print(sed.grid.highlimits)
## ['75840', '72000']
print(sed.boundingboxes)
## [{'nativeSrs': 'EPSG:4326', 'bbox': (-36.0000001152, 14.998958381330127, 43.00104180426987, 90.00000028800001)}]
cx,cy = map(int,sed.grid.highlimits)
bbox = sed.boundingboxes[0]['bbox']
lx,ly,hx,hy = map(float,bbox)
resx,resy = (hx-lx)/cx,(hy-ly)/cy
width = cx/1000
height = cy/1000

We download the data using WCS and save it to a GEOTIFF file.

# get the data
gc = wcs.getCoverage(identifier=bathyml, bbox = requestbbox, coverage=sed, format='GeoTIFF', crs=sed.boundingboxes[0]['nativeSrs'],resx=resx,resy=resy)

# write to a file
fn = clipfile
f = open(fn,'wb')
f.write(gc.read())
## 27648645
f.close()

The data can now be loaded to make a nice bathymetry map using the basemap library.

# load the geotiff file
ds = gdal.Open(clipfile)

# get the dimentions of column and row
nc   = ds.RasterXSize
nr   = ds.RasterYSize

# read elevation data
bathy = ds.ReadAsArray()

# only get positive depth values
bathy[bathy < 0] = 0

# get Longitude and Latitude info
geotransform = ds.GetGeoTransform()
xOrigin      = geotransform[0]
yOrigin      = geotransform[3]
pixelWidth   = geotransform[1]
pixelHeight  = geotransform[5]

# enerate Longitude and Latitude array
lons = xOrigin + np.arange(0, nc)*pixelWidth
lats = yOrigin + np.arange(0, nr)*pixelHeight

# Plot setup
fig= plt.figure(figsize=(10,10))
ax = plt.subplot(111,aspect = 'equal')
plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0, hspace=0)

#Map setup
map = Basemap(resolution='f', projection='cyl', llcrnrlon=requestbbox[0],llcrnrlat=requestbbox[1], urcrnrlon=requestbbox[2], urcrnrlat=requestbbox[3]) #EDIT
parallels = np.arange(-90,90,0.5)
meridians = np.arange(0,360,0.5)
map.drawparallels(parallels,labels=[1,0,0,0],color='w', fontsize=10, fontweight='bold')
meri = map.drawmeridians(meridians,labels=[0,0,0,1],color='w', fontsize=10, fontweight='bold')

#Load colormap and setup elevation contour levels
cmap=plt.get_cmap("GnBu")
keys = np.linspace(round(bathy.min(),-1), round(bathy.max(),-1), 36)

#Contour plot
x, y = map(*np.meshgrid(lons, lats))
cs=map.contourf(x, y, bathy, keys, cmap=cmap)

#add land area
map.fillcontinents(color="#FFDDCC", lake_color='#DDEEFF')
map.drawmapboundary(fill_color="#DDEEFF")
map.drawcoastlines()

map.drawparallels(parallels,labels=[1,0,0,0],color='k', fontsize=10, fontweight='bold')
meri = map.drawmeridians(meridians,labels=[0,0,0,1],color='k', fontsize=10, fontweight='bold')

cb = map.colorbar(cs, 'bottom', size='5%', pad='10%')

cb.set_label('Elevation (m)', fontsize=12, fontweight='bold')
cb.ax.tick_params(labelsize=10)

plt.show()