Versions Compared


  • This line was added.
  • This line was removed.
  • Formatting was changed.


  1. go to the catalogue of EMODnet to find the data sets, these are:
    1. bathymetry
    2. biotic observations (without some specific knowledge it is hard to find data set(s) where abundance information of Amphiura filiformis can be found)
  2. Go to the specific sub portals via the links presented above and find out which services are available, these are:
    1. for bathymetry the following services are made available
      WCS gives you the real data. The example will be worked out for this service. Find out more about WCS at the page WCS primer.
    2. for biotic observations the challenge is to find abundance information on the VLIZ geoserver
      Geoserver can advertise WMS, WFS and WCS information, in case of species distribution in relation to depth the real observations are the way to go. Read more about WFS at the WFS primer page. Checking the complete geoserver and filtering on Eurobis it reveals the best data set to be used is Aggregated points (stations) of EurOBIS occurences.
  3. At this point insight is gained in the available sources and possible data sets that can be used. Close investigation reveals however that the links offered for Bathymetry are not really usable. Reverse engineering via the download area of interest button on the EMODnet bathymetry portal reveals that the data set with the real data is not advertised. This is shown via the get capabilities request of the WCS services provided.

    Code Block
    titleWCS GetCapabilities

    this gives:

    The developer tools of the browser used on the bathymetry data portal gives information about the layer used while downloading an area of interest. The developer tools point towards the coverage emodnet:mean which however is not mentioned in the getCapabilities. This is a bit confusing, especially if you go and expect that the OWS library (used below) results in an error if you do it the way that is advertised according to the OWS documentation (which is not very helpfull).

  4. Now a list of exact data sources is know, these are:
    1. for bathymetry ( coverage emodnet:mean)
    2. for biota ( layername=Eurobis:eurobis_points)
  5. On the primer pages present code snippets for python. These code snippets are adjusted to the particular example presented here.

    Code Block
    titleWCS code snippet
    # import packages
    import os
    from owslib.wcs import WebCoverageService
    from osgeo import gdal
    # define the connection	
    url = ''
    wcs = WebCoverageService(url,version='1.0.0',timeout=320)
    # define variables
    clipfile = r'..\temp.tif'
    requestbbox = (2.097,52.715,4.277,53.935)
    layer = 'emodnet:mean_atlas_land'
    # get the data
    bathyml = 'emodnet:mean'
    sed = wcs[layer] #this is necessary to get essential metadata from the layers advertised
    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
    gc = wcs.getCoverage(identifier=bathyml,
    fn = clipfile
    f = open(fn,'wb')
    filetiff = gdal.Open(clipfile)
    theImage = filetiff.GetRasterBand(1).ReadAsArray()
  6. for querying the WFS following code is used

    Code Block
    titleWFS querying
    # import the packages
    import pandas
    import sys
    import os
    from owslib.fes import *
    from owslib.etree import etree
    from owslib.wfs import WebFeatureService
    # define the input variables
    requestbbox = (2.097,52.715,4.277,53.935)
    fn = r'..\afile.csv'
    wfs = ''  
    layer = 'Eurobis:eurobis_points'
    afilter = PropertyIsEqualTo(propertyname='Eurobis:ScientificName', literal='Amphiura filiformis')
    # define the link to the service and carry out the request
    wfs11 = WebFeatureService(url=wfs, version='1.1.0',timeout=320)
    response = wfs11.getfeature(typename=layer, bbox=bbxrequestbbox,srsname='urn:x-ogc:def:crs:EPSG:4326',outputFormat='csv')   
    if os.path.isfile(fn):
    out = open(fn, 'wb')
    # define pandas dataframe for the output
    df = pandas.read_csv(fn,header=0,sep=',')
  7. Now the result (if there is data in the boundingbox, there is no testing involved, thus the result can be empty, this is not tested) is there, in 2 objects (Pandas DataFrame with observation information and a Numpy array with bathymetry). This can be plotted in a 3D Scatter plot using the Matplotlib. See next code.

    Code Block
    titlePlotting the data
    img = theImage (see above)
    bbox = requestbbox (see above)
    def plotraster(img,df,resx,resy,bbox):
        from mpl_toolkits.mplot3d import Axes3D
        import matplotlib.pyplot as plt
        from matplotlib import cm
        from matplotlib.ticker import LinearLocator, FormatStrFormatter
        import numpy as np
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        #plot the vectors
        xv = df['Longitude']
        yv = df['Latitude']
        zv = df['MaximumDepth']
        # Make data.
        X = np.linspace(bbox[0],bbox[2], img.shape[1])
        Y = np.linspace(bbox[1],bbox[3], img.shape[0])
        X, Y = np.meshgrid(X, Y)
        #R = np.sqrt(X**2 + Y**2)
        #Z = np.sin(R)
        Z = img
        # Plot the surface.
        surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,
                               linewidth=0, antialiased=False,alpha=0.5)
        # Customize the z axis.
        ax.set_zlim(Z.min().round(), Z.max().round())
        # Add a color bar which maps values to colors.
        fig.colorbar(surf, shrink=0.5, aspect=5)
  8. The result is visualised in the following graph