Date: Thu, 28 Mar 2024 15:35:20 +0100 (CET) Message-ID: <1399130194.2564.1711636520386@v-public003.directory.intra> Subject: Exported From Confluence MIME-Version: 1.0 Content-Type: multipart/related; boundary="----=_Part_2563_225353737.1711636520385" ------=_Part_2563_225353737.1711636520385 Content-Type: text/html; charset=UTF-8 Content-Transfer-Encoding: quoted-printable Content-Location: file:///C:/exported.html
At the moment of writing this page, during the 2nd meeting of the Techni=
cal Work Group of EMODnet in the Sea Aquarium of Genova in =
July 2017, discussions about improvements on the services are the main topi=
c. Until now the different EMODnet lots have carried out a great job to col=
lect and harmonize all kinds of data sets. Now it is the time to invest in =
harmonizing all portals so end users can find data a bit better.
At this technical meeting Deltares gave a presentation where a specific us=
e case presented got much attention. The Use Case was "I want to visualise =
the relation between the sea floor depth and the occurrence of a specific a=
nimal (example is for Amphiura filiformis)". Some constraints were=
formulated, these are:
A step wise approach was presented to work out the use case. The steps f= ollowed are:
At this point insight is gained in the available sources and possibl= e 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 po= rtal reveals that the data set with the real data is not advertised. This i= s shown via the get capabilities request of the WCS services provided.
http://o= ws.emodnet-bathymetry.eu/wcs?service=3DWCS&request=3DGetCapabilities&am= p;version=3D1.0.0
this gives:
The developer tools of the browser used on the bathymetry d=
ata portal gives information about the layer used while downloading an =
area of interest. The developer tools point towards the coverage emodnet:me=
an which however is not mentioned in the getCapabilities. This is a bit con=
fusing, especially if you go and expect that the OWS library (used below) r=
esults in an error if you do it the way that is advertised according to the=
OWS documentation (which is not very helpfull).
On the primer pages present code snippets for python. These code sni= ppets are adjusted to the particular example presented here.
# import = packages import os from owslib.wcs import WebCoverageService from osgeo import gdal =20 # define the connection=09 url =3D 'http://ows.emodnet-bathymetry.eu/wcs?' wcs =3D WebCoverageService(url,version=3D'1.0.0',timeout=3D320) wcs.identification.type wcs.identification.title =20 # define variables clipfile =3D r'..\temp.tif' requestbbox =3D (2.097,52.715,4.277,53.935) layer =3D 'emodnet:mean_atlas_land' =20 # get the data bathyml =3D 'emodnet:mean' sed =3D wcs[layer] #this is necessary to get essential metadata from the la= yers advertised sed.keywords sed.grid.highlimits sed.boundingboxes cx, cy =3D map(int,sed.grid.highlimits) bbox =3D sed.boundingboxes[0]['bbox'] lx,ly,hx,hy =3D map(float,bbox) resx, resy =3D (hx-lx)/cx, (hy-ly)/cy width =3D cx/1000 height =3D cy/1000 gc =3D wcs.getCoverage(identifier=3Dbathyml, =09=09 bbox=3Drequestbbox, =09=09 coverage=3Dsed, =09=09 format=3D'GeoTIFF', =09=09 crs=3Dsed.boundingboxes[0]['nativeSrs'], =09=09 width=3Dwidth, =09=09 height=3Dheight) fn =3D clipfile f =3D open(fn,'wb') f.write(gc.read()) f.close() filetiff =3D gdal.Open(clipfile) theImage =3D filetiff.GetRasterBand(1).ReadAsArray() os.unlink(clipfile)
for querying the WFS following code is used
# import = the packages import pandas import sys import os from owslib.fes import * from owslib.etree import etree from owslib.wfs import WebFeatureService =20 # define the input variables requestbbox =3D (2.097,52.715,4.277,53.935) fn =3D r'..\afile.csv' wfs =3D 'http://geo.vliz.be/geoserver/Eurobis/ows?' =20 layer =3D 'Eurobis:eurobis_points' afilter =3D PropertyIsEqualTo(propertyname=3D'Eurobis:ScientificName', lite= ral=3D'Amphiura filiformis') # define the link to the service and carry out the request wfs11 =3D WebFeatureService(url=3Dwfs, version=3D'1.1.0',timeout=3D320) response =3D wfs11.getfeature(typename=3Dlayer, bbox=3Drequestbbox,srsname= =3D'urn:x-ogc:def:crs:EPSG:4326',outputFormat=3D'csv') =20 if os.path.isfile(fn): os.unlink(fn) out =3D open(fn, 'wb') out.write(response.read()) out.close() =20 # define pandas dataframe for the output df =3D pandas.read_csv(fn,header=3D0,sep=3D',')
Now the result (if there is data in the boundingbox, there is no tes= ting involved, thus the result can be empty, this is not tested) is there, = in 2 objects (Pandas DataFrame with observation information and a Numpy arr= ay with bathymetry). This can be plotted in a 3D Scatter plot using the Mat= plotlib. See next code.
img =3D t= heImage (see above) bbox =3D 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 =20 fig =3D plt.figure() ax =3D fig.gca(projection=3D'3d') #plot the vectors xv =3D df['Longitude'] yv =3D df['Latitude'] zv =3D df['MaximumDepth'] ax.scatter(xv,yv,zv,label=3Ddf['ObservedIndividualCount']) =20 # Make data. X =3D np.linspace(bbox[0],bbox[2], img.shape[1]) Y =3D np.linspace(bbox[1],bbox[3], img.shape[0]) X, Y =3D np.meshgrid(X, Y) #R =3D np.sqrt(X**2 + Y**2) #Z =3D np.sin(R) Z =3D img # Plot the surface. surf =3D ax.plot_surface(X, Y, Z, cmap=3Dcm.coolwarm, linewidth=3D0, antialiased=3DFalse,alpha=3D0.5) =20 # Customize the z axis. ax.set_zlim(Z.min().round(), Z.max().round()) ax.zaxis.set_major_locator(LinearLocator(10)) ax.zaxis.set_major_formatter(FormatStrFormatter('%.f')) =20 # Add a color bar which maps values to colors. fig.colorbar(surf, shrink=3D0.5, aspect=3D5) plt.savefig(r'd:\temp\emodnet\test.png',dpi=3D600) =20 plt.show()
The result is visualised in the following graph