Download data from DHI Altimetry API¶

Altimetry data (water level, significant wave height and wind speed) obtained from satellites are very useful for validation of models as the data are available globally since 1985 and up to 12 hour before now.

DHI has an altimetry portal with an api where you can download the data.

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from fmskill.data import DHIAltimetryRepository
from fmskill import TrackObservation
from mikeio import eum
%matplotlib inline

In [2]:

api_key = os.environ["DHI_ALTIMETRY_API_KEY"]
repo = DHIAltimetryRepository(api_key)

In [3]:

repo.get_satellites()

Out[3]:

	long_name
short_name
gs	Geosat
e1	ERS-1
tx	TOPEX
pn	Poseidon
e2	ERS-2
g1	GFO
j1	Jason-1
n1	Envisat
j2	Jason-2
c2	CryoSat-2
sa	SARAL
j3	Jason-3
3a	Sentinel-3A
3b	Sentinel-3B

In [4]:

repo.get_observation_stats()

Out[4]:

	min_date	max_date	count
short_name
c2	2010-07-16 00:12:43.416148	2021-07-01 07:50:33.520384	240692489
n1	2002-05-14 18:32:32.755126	2012-04-08 10:54:23.635591	257633695
j1	2002-01-15 06:07:31.458256	2013-06-21 00:56:32.487362	224260047
sa	2013-03-14 05:39:27.811666	2021-06-30 23:04:13.542166	173100732
j2	2008-07-04 12:19:19.570865	2019-10-01 08:04:57.954027	225297275
pn	1992-10-01 16:45:03.540531	2002-07-12 14:28:19.397424	15718916
gs	1985-03-31 00:00:00.861264	1989-12-30 14:57:52.901281	75115692
j3	2016-02-12 01:11:09.095919	2021-07-01 04:30:25.807815	114658488
e1	1991-08-01 00:25:04.644073	1996-06-02 21:59:26.618659	83641070
tx	1992-09-25 05:17:45.663910	2005-10-08 23:49:39.862525	262462616
e2	1995-04-29 09:39:00.086310	2011-07-04 11:28:00.141221	168846224
3b	2018-05-08 07:42:54.000000	2021-07-01 06:13:07.000000	66721425
3a	2016-03-01 09:20:18.000000	2021-07-01 06:53:02.000000	113451524
g1	2000-01-07 17:28:19.663586	2008-09-17 20:30:05.226965	150233671

In [5]:

repo.plot_observation_stats();

In [6]:

repo.time_of_newest_data

Out[6]:

Timestamp('2021-07-01 07:50:33.520384')

Temporal coverage¶

For a specific area and temporal selection, an overview of the temporal coverage can be obtained with the get_daily_count() method.

In [7]:

df = repo.get_daily_count("lon=10.9&lat=55.9&radius=10.0", start_time="2021")
df.plot(marker="*", figsize=(10,6));

Spatial coverage¶

An overview of the spatial data coverage for a specified area and time period can be obtained with the get_spatial_coverage() method

In [8]:

area = "bbox=10.0,54.0,12.0,58.0"
gdf = repo.get_spatial_coverage(area, start_time="2021")

In [9]:

gdf.plot('count');

In [10]:

# Download coastline and countries from public source
import geopandas
coastline = geopandas.read_file("https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_coastline.geojson")
countries = geopandas.read_file("https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_admin_0_countries.geojson")

In [11]:

ax = gdf.plot('count', legend=True, figsize=(8,8), legend_kwds={'label': "# obs", 'orientation': "horizontal"}, cmap='YlOrRd')
countries.plot(ax=ax, color='0.8')
ax.set_xlim((8,13))
ax.set_ylim((54,58.3))
ax.set_aspect(1.0 / np.cos(np.pi * 56 / 180))

Download the data¶

The actual data is downloaded using the get_altimetry_data() method. Let's download a month's data a circle with 100km radius in the North Sea.

In [12]:

data = repo.get_altimetry_data(area="lon=2.9&lat=55.9&radius=100", start_time="2019-10-1", 
end_time="2019-11-1")
data.df.head()

Succesfully retrieved 1006 records from API in 0.83 seconds

Out[12]:

	lon	lat	water_level	adt	wind_speed	wind_speed_abdalla	wind_speed_abdalla_adjusted	swh	swh_adjusted	distance_from_land	water_depth	satellite	quality	wl_rms	swh_rms
date
2019-10-01 01:53:30.689	3.780429	55.203699	0.0077	0.0597	5.16	5.834183	5.783530	1.975	1.942423	221652.719762	-41.000000	c2	0	0.046	0.256
2019-10-01 01:53:31.633	3.770297	55.260700	0.0287	0.0654	5.21	5.886469	5.836653	1.922	1.892914	227674.772499	-33.931813	c2	0	0.048	0.333
2019-10-01 01:53:32.576	3.760148	55.317701	0.0558	0.0801	5.40	6.101864	6.055494	1.864	1.838815	233715.397905	-28.000000	c2	0	0.038	0.273
2019-10-01 01:53:33.519	3.749980	55.374701	0.0496	0.0645	5.56	6.298769	6.255549	1.962	1.930273	239773.215489	-28.000000	c2	0	0.050	0.374
2019-10-01 01:53:34.463	3.739795	55.431700	0.0752	0.0815	5.33	6.019907	5.972225	1.986	1.952707	245847.148937	-31.000000	c2	0	0.048	0.311

In [13]:

data.query_params

Out[13]:

{'lon': '2.9',
 'lat': '55.9',
 'radius': '100',
 'start_date': '20191001',
 'end_date': '20191101'}

Quality flag¶

Each data point has a quality flag:

Good
Acceptable
Bad

Let's keep only the good data.

In [14]:

data.df = data.df[data.df.quality==0]

In [15]:

data.df[['water_level']].plot(style="*");

In [16]:

data.plot_map()
ax = plt.gca()
countries.plot(ax=ax, color='0.8')
ax.set_xlim((-1,10))
ax.set_ylim((52,60))
ax.set_aspect(1.0 / np.cos(np.pi * 56 / 180))

fmskill.TrackObservation object¶

In [17]:

obs = TrackObservation(data.df, item=7, name='Altimetry c2, 3a, sa, j3')
obs.itemInfo = eum.ItemInfo(eum.EUMType.Significant_wave_height)
obs

Out[17]:

TrackObservation: Altimetry c2, 3a, sa, j3, n=715

In [18]:

obs.hist();

It is now ready for comparison with you model results...

Get all observations in model domain¶

In [19]:

from fmskill import ModelResult
fn = '../tests/testdata/SW/HKZN_local_2017_DutchCoast.dfsu'
mr = ModelResult(fn, name='MIKE21SW', item=0)
mr.dfs

Out[19]:

Dfsu2D
Number of elements: 958
Number of nodes: 570
Projection: LONG/LAT
Number of items: 15
Time: 23 steps with dt=10800.0s
      2017-10-27 00:00:00 -- 2017-10-29 18:00:00

In [20]:

# string representation of domain boundary polyline
xy = mr.dfs.boundary_polylines.exteriors[0].xy
SW_domain=np.array2string(xy.flatten(), formatter={'float_kind':lambda x: "%.4f" % x}, max_line_width=1000000, separator=',').replace('[','polygon=').replace(']','')

In [21]:

start = mr.start_time
end = mr.end_time
data = repo.get_altimetry_data(area=SW_domain, start_time=start, end_time=end)
data.df.head()

Succesfully retrieved 367 records from API in 0.49 seconds

Out[21]:

	lon	lat	water_level	adt	wind_speed	wind_speed_abdalla	wind_speed_abdalla_adjusted	swh	swh_adjusted	distance_from_land	water_depth	satellite	quality	wl_rms	swh_rms
date
2017-10-27 10:45:19	1.262463	55.298639	0.3778	0.1940	9.23	9.581609	9.600170	2.255	2.240740	148856.276678	-67.506027	3a	0	0.0587	0.400
2017-10-27 10:45:20	1.230990	55.241321	0.4375	0.2411	9.36	9.723277	9.743679	2.130	2.122240	143450.079792	-64.900677	3a	0	0.0424	0.464
2017-10-27 10:45:21	1.199594	55.183996	0.4489	0.2399	9.19	9.581609	9.600170	2.279	2.263492	138052.640961	-63.659107	3a	0	0.0455	0.349
2017-10-27 10:45:22	1.168277	55.126665	0.4711	0.2490	9.17	9.546213	9.564314	2.129	2.121292	132754.802194	-64.202776	3a	0	0.0413	0.286
2017-10-27 10:45:23	1.137037	55.069327	0.4678	0.2323	9.09	9.440084	9.456805	2.084	2.078632	127587.167133	-58.527079	3a	0	0.0600	0.361

In [22]:

data.satellites

Out[22]:

['3a', 'c2', 'j2', 'j3', 'sa']

In [23]:

data.df.swh.plot(style='o');

In [24]:

ax = mr.dfs.plot(plot_type='outline_only', figsize=(10,8))
countries.plot(ax=ax, color='0.8')
df = data.df
df[df.quality==0].plot.scatter(x='lon', y='lat', c='swh', cmap='YlOrRd', ax=ax);
df[df.quality==1].plot.scatter(x='lon', y='lat', c='blue', ax=ax);
df[df.quality==2].plot.scatter(x='lon', y='lat', c='black', ax=ax);

In [25]:

data.print_records_per_satellite()

For the selected area between 2017-10-27 10:45:19 and 2017-10-28 20:36:40.787000:
Satellite 3a has 64 records
Satellite c2 has 83 records
Satellite j2 has 87 records
Satellite j3 has 108 records
Satellite sa has 25 records

In [26]:

data.to_dfs0('alti_data.dfs0')

Clean up¶

In [27]:

import os
os.remove('alti_data.dfs0')