import mikeio
+
-
ds = mikeio.read("../data/wind_north_sea.dfsu", items="Wind speed")
ds
@@ -398,7 +398,7 @@
Dfsu - 2D interpolation
0: Wind speed <Wind speed> (meter per sec)
+
da = ds.Wind_speed
da.plot();
@@ -416,7 +416,7 @@ The interpolated data is then saved to a dfs2 file.
-
Interpolate to gridThen interpolate all data to the new grid and plot.
+
-
g = da.geometry.get_overset_grid(dx=0.1)
g
@@ -426,7 +426,7 @@
Interpolate to grid
+
-
da_grid = da.interp_like(g)
da_grid
@@ -437,7 +437,7 @@
Interpolate to grid
+
da_grid.plot();
@@ -450,14 +450,14 @@ Interpolate to grid
Interpolate to grid
Save to dfs2 file
-
+
da_grid.to_dfs("wind_north_sea_interpolated.dfs2")
Save to NetCDF
-
+
xr_da = da_grid.to_xarray()
xr_da.to_netcdf("wind_north_sea_interpolated.nc")
@@ -478,7 +478,7 @@ Save to GeoTiff
This section requires the rasterio package.
-
+
import numpy as np
import rasterio
from rasterio.transform import from_origin
@@ -502,7 +502,7 @@ Save to GeoTiff
Interpolate to other mesh
Interpolate the data from this coarse mesh onto a finer resolution mesh
-
+
msh = mikeio.Mesh('../data/north_sea_2.mesh')
msh
@@ -512,7 +512,7 @@ Interpolate to other mesh
projection: LONG/LAT
-
+
dsi = da.interp_like(msh)
dsi
@@ -523,7 +523,7 @@ Interpolate to other mesh
geometry: Dfsu2D (2259 elements, 1296 nodes)
-
+
da[0].plot(figsize=(9,7), show_mesh=True);
@@ -533,7 +533,7 @@ Interpolate to other mesh
-
+
dsi[0].plot(figsize=(9,7), show_mesh=True);
@@ -545,14 +545,14 @@ Interpolate to other mesh
Note: 3 of the new elements are outside the original mesh and data are therefore NaN by default
-
+
nan_elements = np.where(np.isnan(dsi[0].to_numpy()))[0]
nan_elements
array([ 249, 451, 1546])
-
+
da.geometry.contains(msh.element_coordinates[nan_elements,:2])
array([False, False, False])
@@ -561,10 +561,10 @@
We can force extrapolation to avoid the NaN values
-
+
dat_interp = da.interp_like(msh, extrapolate=True)
-
+
n_nan_elements = np.sum(np.isnan(dat_interp.values))
n_nan_elements
@@ -577,14 +577,14 @@ Interpola
We want to interpolate scatter data onto an existing mesh and create a new dfsu with the interpolated data.
This uses lower level private utility methods not part of the public API.
Interpolating from scatter data will soon be possible in a simpler way.
-
+
from mikeio.spatial._utils import dist_in_meters
from mikeio._interpolation import get_idw_interpolant
-
+
dfs = mikeio.open("../data/wind_north_sea.dfsu")
-
+
dfs.geometry.plot.mesh();
@@ -594,7 +594,7 @@ Interpola
-
+
# scatter data: x,y,value for 4 points
scatter= np.array([[1,50,1], [4, 52, 3], [8, 55, 2], [-1, 55, 1.5]])
scatter
@@ -611,35 +611,35 @@ Interpola
calc IDW interpolatant weights
Interpolate
-
+
dist = dist_in_meters(scatter[:,:2], dfs.element_coordinates[0,:2])
dist
array([4.00139539, 3.18881018, 6.58769411, 2.69722991])
-
+
w = get_idw_interpolant(dist, p=2)
w
array([0.19438779, 0.30607974, 0.07171749, 0.42781498])
-
+
np.dot(scatter[:,2], w) # interpolated value in element 0
1.8977844597276883
Let’s do the same for all points in the mesh and plot in the end
-
+
dati = np.zeros((1,dfs.n_elements))
for j in range(dfs.n_elements):
dist = dist_in_meters(scatter[:,:2], dfs.element_coordinates[j,:2])
w = get_idw_interpolant(dist, p=2)
dati[0,j] = np.dot(scatter[:,2], w)
-
+
da = mikeio.DataArray(data=dati, geometry=dfs.geometry, time=dfs.start_time)
da
@@ -650,7 +650,7 @@ Interpola
geometry: Dfsu2D (958 elements, 570 nodes)
-
+
da.plot(title="Interpolated scatter data");
@@ -660,13 +660,13 @@ Interpola
-
+
da.to_dfs("interpolated_scatter.dfsu")
Clean up
-
+
import os
os.remove("wind_north_sea_interpolated.dfs2")
diff --git a/examples/Generic.html b/examples/Generic.html
index 32a1b96ce..838eb59d0 100644
--- a/examples/Generic.html
+++ b/examples/Generic.html
@@ -265,7 +265,7 @@ Generic dfs processing
@@ -394,7 +394,7 @@ Generic dfs processing
quantile: Create temporal quantiles of dfs file
-
+
import matplotlib.pyplot as plt
import mikeio
import mikeio.generic
@@ -402,7 +402,7 @@ Generic dfs processing
Concatenation
Take a look at these two files with overlapping timesteps.
-
+
t1 = mikeio.read("../data/tide1.dfs1")
t1
@@ -414,7 +414,7 @@ Concatenation
0: Level <Water Level> (meter)
-
+
t2 = mikeio.read("../data/tide2.dfs1")
t2
@@ -427,7 +427,7 @@ Concatenation
Plot one of the points along the line.
-
+
plt.plot(t1.time,t1[0].isel(x=1).values, label="File 1")
plt.plot(t2.time,t2[0].isel(x=1).values,'k+', label="File 2")
plt.legend()
@@ -439,15 +439,15 @@ Concatenation
-
+
mikeio.generic.concat(infilenames=["../data/tide1.dfs1",
"../data/tide2.dfs1"],
outfilename="concat.dfs1")
- 0%| | 0/2 [00:00<?, ?it/s]100%|██████████| 2/2 [00:00<00:00, 627.04it/s]
+ 0%| | 0/2 [00:00<?, ?it/s]100%|██████████| 2/2 [00:00<00:00, 525.08it/s]
-
+
c = mikeio.read("concat.dfs1")
c[0].isel(x=1).plot()
c
@@ -471,16 +471,16 @@ Concatenation
Difference between two files
Take difference between two dfs files with same structure - e.g. to see the difference in result between two calibration runs
-
+
fn1 = "../data/oresundHD_run1.dfsu"
fn2 = "../data/oresundHD_run2.dfsu"
fn_diff = "oresundHD_difference.dfsu"
mikeio.generic.diff(fn1, fn2, fn_diff)
- 0%| | 0/5 [00:00<?, ?it/s]100%|██████████| 5/5 [00:00<00:00, 2690.73it/s]
+ 0%| | 0/5 [00:00<?, ?it/s]100%|██████████| 5/5 [00:00<00:00, 2250.16it/s]
-
+
_, ax = plt.subplots(1,3, sharey=True, figsize=(12,5))
da = mikeio.read(fn1, time=-1)[0]
da.plot(vmin=0.06, vmax=0.27, ax=ax[0], title='run 1')
@@ -504,11 +504,11 @@ Extract time s
time slice by specifying start and/or end
specific items
-
+
infile = "../data/tide1.dfs1"
mikeio.generic.extract(infile, "extracted.dfs1", start='2019-01-02')
-
+
e = mikeio.read("extracted.dfs1")
e
@@ -520,11 +520,11 @@ Extract time s
0: Level <Water Level> (meter)
-
+
infile = "../data/oresund_vertical_slice.dfsu"
mikeio.generic.extract(infile, "extracted.dfsu", items='Salinity', end=-2)
-
+
e = mikeio.read("extracted.dfsu")
e
@@ -545,7 +545,7 @@ Extract time s
Scaling
Adding a constant e.g to adjust datum
-
+
ds = mikeio.read("../data/gebco_sound.dfs2")
ds.Elevation[0].plot();
@@ -556,23 +556,23 @@ Scaling
-
+
ds['Elevation'][0,104,131].to_numpy()
-1.0
This is the processing step.
-
+
mikeio.generic.scale("../data/gebco_sound.dfs2",
"gebco_sound_local_datum.dfs2",
offset=-2.1
)
- 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1047.27it/s]
+ 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1284.23it/s]
-
+
ds2 = mikeio.read("gebco_sound_local_datum.dfs2")
ds2['Elevation'][0].plot()
@@ -583,7 +583,7 @@ Scaling
-
+
ds2['Elevation'][0,104,131].to_numpy()
-3.1
@@ -591,7 +591,7 @@ Scaling
Spatially varying correction
-
+
import numpy as np
factor = np.ones_like(ds['Elevation'][0].to_numpy())
factor.shape
@@ -600,7 +600,7 @@ Spatially var
Add some spatially varying factors, exaggerated values for educational purpose.
-
+
factor[:,0:100] = 5.3
factor[0:40,] = 0.1
factor[150:,150:] = 10.7
@@ -615,7 +615,7 @@ Spatially var
The 2d array must first be flipped upside down and then converted to a 1d vector using numpy.ndarray.flatten to match how data is stored in dfs files.
-
+
factor_ud = np.flipud(factor)
factor_vec = factor_ud.flatten()
mikeio.generic.scale("../data/gebco_sound.dfs2",
@@ -623,10 +623,10 @@ Spatially var
factor=factor_vec
)
- 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1223.54it/s]
+ 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1265.25it/s]
-
+
ds3 = mikeio.read("gebco_sound_spatial.dfs2")
ds3.Elevation[0].plot();
@@ -641,15 +641,15 @@ Spatially var
Time average
-
+
fn = "../data/NorthSea_HD_and_windspeed.dfsu"
fn_avg = "Avg_NorthSea_HD_and_windspeed.dfsu"
mikeio.generic.avg_time(fn, fn_avg)
- 0%| | 0/66 [00:00<?, ?it/s]100%|██████████| 66/66 [00:00<00:00, 17965.09it/s]
+ 0%| | 0/66 [00:00<?, ?it/s]100%|██████████| 66/66 [00:00<00:00, 15917.66it/s]
-
+
ds = mikeio.read(fn)
ds.mean(axis=0).describe() # alternative way of getting the time average
@@ -713,7 +713,7 @@ Time average
-
+
ds_avg = mikeio.read(fn_avg)
ds_avg.describe()
@@ -781,12 +781,12 @@ Time average
Quantile
Example that calculates the 25%, 50% and 75% percentile for all items in a dfsu file.
-
+
fn = "../data/NorthSea_HD_and_windspeed.dfsu"
fn_q = "Q_NorthSea_HD_and_windspeed.dfsu"
mikeio.generic.quantile(fn, fn_q, q=[0.25,0.5,0.75])
-
+
ds = mikeio.read(fn_q)
ds
@@ -803,7 +803,7 @@ Quantile
5: Quantile 0.75, Wind speed <Wind speed> (meter per sec)
-
+
da_q75 = ds["Quantile 0.75, Wind speed"]
da_q75.plot(title="75th percentile, wind speed", label="m/s")
@@ -817,7 +817,7 @@ Quantile
Clean up
-
+
import os
os.remove("concat.dfs1")
os.remove("oresundHD_difference.dfsu")
diff --git a/examples/Time-interpolation.html b/examples/Time-interpolation.html
index 420005818..d1ffd0a12 100644
--- a/examples/Time-interpolation.html
+++ b/examples/Time-interpolation.html
@@ -263,7 +263,7 @@ Time interpolation
@@ -376,11 +376,11 @@ Time interpolation
-
+
import numpy as np
import mikeio
-
+
ds = mikeio.read("../data/waves.dfs2")
ds
@@ -397,7 +397,7 @@ Time interpolation
Interpolate to specific timestep
A common use case is to interpolate to a shorter timestep, in this case 1h.
-
+
ds_h = ds.interp_time(3600)
ds_h
@@ -412,14 +412,14 @@ Interpola
And to store the interpolated data in a new file.
-
+
ds_h.to_dfs("waves_3h.dfs2")
Interpolate to time axis of another dataset
Read some non-equidistant data typically found in observed data.
-
+
ts = mikeio.read("../data/waves.dfs0")
ts
@@ -434,10 +434,10 @@
+
dsi = ds.interp_time(ts)
-
+
dsi.time
DatetimeIndex(['2004-01-01 01:00:00', '2004-01-01 02:00:00',
@@ -455,13 +455,13 @@
+
dsi["Sign. Wave Height"].shape
(24, 31, 31)
-
+
ax = dsi["Sign. Wave Height"].sel(x=250, y=1200).plot(marker='+')
ts["Sign. Wave Height"].plot(ax=ax,marker='+')
@@ -479,7 +479,7 @@ Model validation
In the example below we calculate this metric using the model data interpolated to the observed times.
For a more elaborate model validation library which takes care of these things for you as well as calculating a number of relevant metrics, take a look at `ModelSkill.
Use np.nanmean to skip NaN.
-
+
ts["Sign. Wave Height"]
<mikeio.DataArray>
@@ -490,7 +490,7 @@ Model validation
values: [0.06521, 0.06771, ..., 0.0576]
-
+
dsi["Sign. Wave Height"].sel(x=250, y=1200)
<mikeio.DataArray>
@@ -501,7 +501,7 @@ Model validation
values: [0.0387, 0.03939, ..., nan]
-
+
diff = (ts["Sign. Wave Height"] - dsi["Sign. Wave Height"].sel(x=250, y=1200))
diff.plot()
@@ -512,7 +512,7 @@ Model validation
-
+
mae = np.abs(diff).nanmean().to_numpy()
mae
@@ -522,7 +522,7 @@ Model validation
Clean up
-
+
import os
os.remove("waves_3h.dfs2")
diff --git a/examples/dfs2/bathy.html b/examples/dfs2/bathy.html
index 5960c44c9..a1f61fc0f 100644
--- a/examples/dfs2/bathy.html
+++ b/examples/dfs2/bathy.html
@@ -266,7 +266,7 @@ Dfs2 - Bathymetric data
@@ -377,11 +377,11 @@ Dfs2 - Bathymetric data
GEBCO Compilation Group (2020) GEBCO 2020 Grid (doi:10.5285/a29c5465-b138-234d-e053-6c86abc040b9)
-
+
import xarray
import mikeio
-
+
ds = xarray.open_dataset("../../data/gebco_2020_n56.3_s55.2_w12.2_e13.1.nc")
ds
@@ -753,7 +753,7 @@ Dfs2 - Bathymetric data
Coordinates: (2)
Data variables:
elevation (lat, lon) int16 114kB ...
-Attributes: (8)- Conventions :
- CF-1.6
- title :
- The GEBCO_2020 Grid - a continuous terrain model for oceans and land at 15 arc-second intervals
- institution :
- On behalf of the General Bathymetric Chart of the Oceans (GEBCO), the data are held at the British Oceanographic Data Centre (BODC).
- source :
- The GEBCO_2020 Grid is the latest global bathymetric product released by the General Bathymetric Chart of the Oceans (GEBCO) and has been developed through the Nippon Foundation-GEBCO Seabed 2030 Project. This is a collaborative project between the Nippon Foundation of Japan and GEBCO. The Seabed 2030 Project aims to bring together all available bathymetric data to produce the definitive map of the world ocean floor and make it available to all.
- history :
- Information on the development of the data set and the source data sets included in the grid can be found in the data set documentation available from https://www.gebco.net
- references :
- DOI: 10.5285/a29c5465-b138-234d-e053-6c86abc040b9
- comment :
- The data in the GEBCO_2020 Grid should not be used for navigation or any purpose relating to safety at sea.
- node_offset :
- 1.0
ds.elevation.plot();
@@ -784,7 +784,7 @@ Dfs2 - Bathymetric data
-
+
ds.elevation.sel(lon=12.74792, lat=55.865, method="nearest")
+Attributes: (7)
Check ordering of dimensions, should be (y,x)
-
+
ds.elevation.dims
('lat', 'lon')
-
+
el = ds.elevation.values
el.shape
@@ -1171,37 +1171,37 @@ Dfs2 - Bathymetric data
Check that axes are increasing, S->N W->E
-
+
ds.lat.values[0],ds.lat.values[-1]
(55.20208333333332, 56.29791666666665)
-
+
ds.lat.values[0] < ds.lat.values[-1]
True
-
+
ds.lon.values[0],ds.lon.values[-1]
(12.20208333333332, 13.097916666666663)
-
+
el[0,0] # Bottom left
-8
-
+
el[-1,0] # Top Left
-31
-
+
geometry = mikeio.Grid2D(x=ds.lon.values, y=ds.lat.values, projection="LONG/LAT")
geometry
@@ -1211,7 +1211,7 @@ Dfs2 - Bathymetric data
projection: LONG/LAT
-
+
da = mikeio.DataArray(data=el,
item=mikeio.ItemInfo("Elevation", mikeio.EUMType.Total_Water_Depth),
geometry=geometry,
@@ -1226,7 +1226,7 @@ Dfs2 - Bathymetric data
geometry: Grid2D (ny=264, nx=216)
-
+
da.plot();
@@ -1236,7 +1236,7 @@ Dfs2 - Bathymetric data
-
+
da.plot(cmap='coolwarm', vmin=-100, vmax=100);
@@ -1246,10 +1246,10 @@ Dfs2 - Bathymetric data
-
+
da.to_dfs("gebco.dfs2")
-
+
ds = mikeio.read("gebco.dfs2")
ds.Elevation.plot()
@@ -1262,7 +1262,7 @@ Dfs2 - Bathymetric data
Clean up
-
+
import os
os.remove("gebco.dfs2")
diff --git a/examples/dfs2/gfs.html b/examples/dfs2/gfs.html
index 3c0467a6e..461e0bad8 100644
--- a/examples/dfs2/gfs.html
+++ b/examples/dfs2/gfs.html
@@ -266,7 +266,7 @@ Dfs2 - Meteo data
@@ -380,13 +380,13 @@ Dfs2 - Meteo data
-
+
import xarray
import pandas as pd
import mikeio
The file gfs_wind.nc contains a small sample of the GFS forecast data downloaded via their OpenDAP service
-
+
ds = xarray.open_dataset('../../data/gfs_wind.nc')
ds
@@ -760,30 +760,30 @@ Dfs2 - Meteo data
msletmsl (time, lat, lon) float32 10kB ...
ugrd10m (time, lat, lon) float32 10kB ...
vgrd10m (time, lat, lon) float32 10kB ...
-Attributes: (4)
Running a Mike 21 HD model, needs at least three variables of meteorological forcing * Mean Sea Level Pressure * U 10m * V 10m
Let’s take a look the U 10m
-
+
ds.ugrd10m.isel(time=0).plot();
@@ -797,7 +797,7 @@ Dfs2 - Meteo data
Convert to dfs2
Time
-
+
time = pd.DatetimeIndex(ds.time)
time
@@ -809,36 +809,36 @@ Time
Variable types
-
+
mikeio.EUMType.Air_Pressure
Air Pressure
-
+
mikeio.EUMType.Air_Pressure.units
[hectopascal, millibar]
-
+
mikeio.EUMType.Wind_Velocity
Wind Velocity
-
+
mikeio.EUMType.Wind_Velocity.units
[meter per sec, feet per sec, miles per hour, km per hour, knot]
-
+
mslp = ds.msletmsl.values / 100 # conversion from Pa to hPa
u = ds.ugrd10m.values
v = ds.vgrd10m.values
-
+
geometry = mikeio.Grid2D(x=ds.lon.values, y=ds.lat.values, projection="LONG/LAT")
geometry
@@ -848,14 +848,14 @@ Variable types
projection: LONG/LAT
-
+
from mikeio import ItemInfo, EUMType, EUMUnit
mslp_da = mikeio.DataArray(data=mslp,time=time, geometry=geometry, item=ItemInfo("Mean Sea Level Pressure", EUMType.Air_Pressure, EUMUnit.hectopascal))
u_da = mikeio.DataArray(data=u,time=time, geometry=geometry, item=ItemInfo("Wind U", EUMType.Wind_Velocity, EUMUnit.meter_per_sec))
v_da = mikeio.DataArray(data=v,time=time, geometry=geometry, item=ItemInfo("Wind V", EUMType.Wind_Velocity, EUMUnit.meter_per_sec))
-
+
mds = mikeio.Dataset([mslp_da, u_da, v_da])
mds
@@ -869,11 +869,11 @@ Variable types
2: Wind V <Wind Velocity> (meter per sec)
-
+
mds.to_dfs("gfs.dfs2")
Clean up
-
+
import os
os.remove("gfs.dfs2")
diff --git a/examples/dfs2/index.html b/examples/dfs2/index.html
index 580c5acd0..c5ff0ca9d 100644
--- a/examples/dfs2/index.html
+++ b/examples/dfs2/index.html
@@ -228,7 +228,7 @@ Dfs2 examples
diff --git a/examples/index.html b/examples/index.html
index b21302007..3cc800c42 100644
--- a/examples/index.html
+++ b/examples/index.html
@@ -288,7 +288,7 @@ Examples
@@ -420,7 +420,7 @@ Examples
-
+
Dfs2 - Bathymetric data
@@ -428,7 +428,7 @@ Examples
Convert GEBCO 2020 NetCDF to dfs2
-
+
Dfs2 - Meteo data
@@ -436,7 +436,7 @@ Examples
Conversion of NetCDF from Global Forecasting System to Dfs2
-
+
Dfs2 examples
@@ -447,7 +447,7 @@ Examples
-
+
Dfsu - 2D interpolation
@@ -455,7 +455,7 @@ Examples
Interpolate dfsu data to a grid, save as dfs2 and geotiff. Interpolate dfsu data to another mesh.
-
+
Generic dfs processing
@@ -463,7 +463,7 @@ Examples
Tools and methods that applies to any type of dfs files.
-
+
Time interpolation
@@ -866,8 +866,8 @@ Examples
Save to dfs2 file
-
+
-da_grid.to_dfs("wind_north_sea_interpolated.dfs2")
Save to NetCDF
-
+
@@ -478,7 +478,7 @@
xr_da = da_grid.to_xarray()
xr_da.to_netcdf("wind_north_sea_interpolated.nc")Save to GeoTiff
This section requires the rasterio package.
+
import numpy as np
import rasterio
from rasterio.transform import from_origin
@@ -502,7 +502,7 @@ Save to GeoTiff
Interpolate to other mesh
Interpolate the data from this coarse mesh onto a finer resolution mesh
-
+
msh = mikeio.Mesh('../data/north_sea_2.mesh')
msh
@@ -512,7 +512,7 @@ Interpolate to other mesh
projection: LONG/LAT
-
+
dsi = da.interp_like(msh)
dsi
@@ -523,7 +523,7 @@ Interpolate to other mesh
geometry: Dfsu2D (2259 elements, 1296 nodes)
-
+
da[0].plot(figsize=(9,7), show_mesh=True);
@@ -533,7 +533,7 @@ Interpolate to other mesh
-
+
dsi[0].plot(figsize=(9,7), show_mesh=True);
@@ -545,14 +545,14 @@ Interpolate to other mesh
Note: 3 of the new elements are outside the original mesh and data are therefore NaN by default
-
+
nan_elements = np.where(np.isnan(dsi[0].to_numpy()))[0]
nan_elements
array([ 249, 451, 1546])
-
+
da.geometry.contains(msh.element_coordinates[nan_elements,:2])
array([False, False, False])
@@ -561,10 +561,10 @@
We can force extrapolation to avoid the NaN values
-
+
dat_interp = da.interp_like(msh, extrapolate=True)
-
+
n_nan_elements = np.sum(np.isnan(dat_interp.values))
n_nan_elements
@@ -577,14 +577,14 @@ Interpola
We want to interpolate scatter data onto an existing mesh and create a new dfsu with the interpolated data.
This uses lower level private utility methods not part of the public API.
Interpolating from scatter data will soon be possible in a simpler way.
-
+
from mikeio.spatial._utils import dist_in_meters
from mikeio._interpolation import get_idw_interpolant
-
+
dfs = mikeio.open("../data/wind_north_sea.dfsu")
-
+
dfs.geometry.plot.mesh();
@@ -594,7 +594,7 @@ Interpola
-
+
# scatter data: x,y,value for 4 points
scatter= np.array([[1,50,1], [4, 52, 3], [8, 55, 2], [-1, 55, 1.5]])
scatter
@@ -611,35 +611,35 @@ Interpola
calc IDW interpolatant weights
Interpolate
-
+
dist = dist_in_meters(scatter[:,:2], dfs.element_coordinates[0,:2])
dist
array([4.00139539, 3.18881018, 6.58769411, 2.69722991])
-
+
w = get_idw_interpolant(dist, p=2)
w
array([0.19438779, 0.30607974, 0.07171749, 0.42781498])
-
+
np.dot(scatter[:,2], w) # interpolated value in element 0
1.8977844597276883
Let’s do the same for all points in the mesh and plot in the end
-
+
dati = np.zeros((1,dfs.n_elements))
for j in range(dfs.n_elements):
dist = dist_in_meters(scatter[:,:2], dfs.element_coordinates[j,:2])
w = get_idw_interpolant(dist, p=2)
dati[0,j] = np.dot(scatter[:,2], w)
-
+
da = mikeio.DataArray(data=dati, geometry=dfs.geometry, time=dfs.start_time)
da
@@ -650,7 +650,7 @@ Interpola
geometry: Dfsu2D (958 elements, 570 nodes)
-
+
da.plot(title="Interpolated scatter data");
@@ -660,13 +660,13 @@ Interpola
-
+
da.to_dfs("interpolated_scatter.dfsu")
Clean up
-
+
import os
os.remove("wind_north_sea_interpolated.dfs2")
diff --git a/examples/Generic.html b/examples/Generic.html
index 32a1b96ce..838eb59d0 100644
--- a/examples/Generic.html
+++ b/examples/Generic.html
@@ -265,7 +265,7 @@ Generic dfs processing
@@ -394,7 +394,7 @@ Generic dfs processing
quantile: Create temporal quantiles of dfs file
-
+
import matplotlib.pyplot as plt
import mikeio
import mikeio.generic
@@ -402,7 +402,7 @@ Generic dfs processing
Concatenation
Take a look at these two files with overlapping timesteps.
-
+
t1 = mikeio.read("../data/tide1.dfs1")
t1
@@ -414,7 +414,7 @@ Concatenation
0: Level <Water Level> (meter)
-
+
t2 = mikeio.read("../data/tide2.dfs1")
t2
@@ -427,7 +427,7 @@ Concatenation
Plot one of the points along the line.
-
+
plt.plot(t1.time,t1[0].isel(x=1).values, label="File 1")
plt.plot(t2.time,t2[0].isel(x=1).values,'k+', label="File 2")
plt.legend()
@@ -439,15 +439,15 @@ Concatenation
-
+
mikeio.generic.concat(infilenames=["../data/tide1.dfs1",
"../data/tide2.dfs1"],
outfilename="concat.dfs1")
- 0%| | 0/2 [00:00<?, ?it/s]100%|██████████| 2/2 [00:00<00:00, 627.04it/s]
+ 0%| | 0/2 [00:00<?, ?it/s]100%|██████████| 2/2 [00:00<00:00, 525.08it/s]
-
+
c = mikeio.read("concat.dfs1")
c[0].isel(x=1).plot()
c
@@ -471,16 +471,16 @@ Concatenation
Difference between two files
Take difference between two dfs files with same structure - e.g. to see the difference in result between two calibration runs
-
+
fn1 = "../data/oresundHD_run1.dfsu"
fn2 = "../data/oresundHD_run2.dfsu"
fn_diff = "oresundHD_difference.dfsu"
mikeio.generic.diff(fn1, fn2, fn_diff)
- 0%| | 0/5 [00:00<?, ?it/s]100%|██████████| 5/5 [00:00<00:00, 2690.73it/s]
+ 0%| | 0/5 [00:00<?, ?it/s]100%|██████████| 5/5 [00:00<00:00, 2250.16it/s]
-
+
_, ax = plt.subplots(1,3, sharey=True, figsize=(12,5))
da = mikeio.read(fn1, time=-1)[0]
da.plot(vmin=0.06, vmax=0.27, ax=ax[0], title='run 1')
@@ -504,11 +504,11 @@ Extract time s
time slice by specifying start and/or end
specific items
-
+
infile = "../data/tide1.dfs1"
mikeio.generic.extract(infile, "extracted.dfs1", start='2019-01-02')
-
+
e = mikeio.read("extracted.dfs1")
e
@@ -520,11 +520,11 @@ Extract time s
0: Level <Water Level> (meter)
-
+
infile = "../data/oresund_vertical_slice.dfsu"
mikeio.generic.extract(infile, "extracted.dfsu", items='Salinity', end=-2)
-
+
e = mikeio.read("extracted.dfsu")
e
@@ -545,7 +545,7 @@ Extract time s
Scaling
Adding a constant e.g to adjust datum
-
+
ds = mikeio.read("../data/gebco_sound.dfs2")
ds.Elevation[0].plot();
@@ -556,23 +556,23 @@ Scaling
-
+
ds['Elevation'][0,104,131].to_numpy()
-1.0
This is the processing step.
-
+
mikeio.generic.scale("../data/gebco_sound.dfs2",
"gebco_sound_local_datum.dfs2",
offset=-2.1
)
- 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1047.27it/s]
+ 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1284.23it/s]
-
+
ds2 = mikeio.read("gebco_sound_local_datum.dfs2")
ds2['Elevation'][0].plot()
@@ -583,7 +583,7 @@ Scaling
-
+
ds2['Elevation'][0,104,131].to_numpy()
-3.1
@@ -591,7 +591,7 @@ Scaling
Spatially varying correction
-
+
import numpy as np
factor = np.ones_like(ds['Elevation'][0].to_numpy())
factor.shape
@@ -600,7 +600,7 @@ Spatially var
Add some spatially varying factors, exaggerated values for educational purpose.
-
+
factor[:,0:100] = 5.3
factor[0:40,] = 0.1
factor[150:,150:] = 10.7
@@ -615,7 +615,7 @@ Spatially var
The 2d array must first be flipped upside down and then converted to a 1d vector using numpy.ndarray.flatten to match how data is stored in dfs files.
-
+
factor_ud = np.flipud(factor)
factor_vec = factor_ud.flatten()
mikeio.generic.scale("../data/gebco_sound.dfs2",
@@ -623,10 +623,10 @@ Spatially var
factor=factor_vec
)
- 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1223.54it/s]
+ 0%| | 0/1 [00:00<?, ?it/s]100%|██████████| 1/1 [00:00<00:00, 1265.25it/s]
-
+
ds3 = mikeio.read("gebco_sound_spatial.dfs2")
ds3.Elevation[0].plot();
@@ -641,15 +641,15 @@ Spatially var
Time average
-
+
fn = "../data/NorthSea_HD_and_windspeed.dfsu"
fn_avg = "Avg_NorthSea_HD_and_windspeed.dfsu"
mikeio.generic.avg_time(fn, fn_avg)
- 0%| | 0/66 [00:00<?, ?it/s]100%|██████████| 66/66 [00:00<00:00, 17965.09it/s]
+ 0%| | 0/66 [00:00<?, ?it/s]100%|██████████| 66/66 [00:00<00:00, 15917.66it/s]
-
+
ds = mikeio.read(fn)
ds.mean(axis=0).describe() # alternative way of getting the time average
@@ -713,7 +713,7 @@ Time average
-
+
ds_avg = mikeio.read(fn_avg)
ds_avg.describe()
@@ -781,12 +781,12 @@ Time average
Quantile
Example that calculates the 25%, 50% and 75% percentile for all items in a dfsu file.
-
+
fn = "../data/NorthSea_HD_and_windspeed.dfsu"
fn_q = "Q_NorthSea_HD_and_windspeed.dfsu"
mikeio.generic.quantile(fn, fn_q, q=[0.25,0.5,0.75])
-
+
ds = mikeio.read(fn_q)
ds
@@ -803,7 +803,7 @@ Quantile
5: Quantile 0.75, Wind speed <Wind speed> (meter per sec)
-
+
da_q75 = ds["Quantile 0.75, Wind speed"]
da_q75.plot(title="75th percentile, wind speed", label="m/s")
@@ -817,7 +817,7 @@ Quantile
Clean up
-
+
import os
os.remove("concat.dfs1")
os.remove("oresundHD_difference.dfsu")
diff --git a/examples/Time-interpolation.html b/examples/Time-interpolation.html
index 420005818..d1ffd0a12 100644
--- a/examples/Time-interpolation.html
+++ b/examples/Time-interpolation.html
@@ -263,7 +263,7 @@ Time interpolation
@@ -376,11 +376,11 @@ Time interpolation
-
+
import numpy as np
import mikeio
-
+
ds = mikeio.read("../data/waves.dfs2")
ds
@@ -397,7 +397,7 @@ Time interpolation
Interpolate to specific timestep
A common use case is to interpolate to a shorter timestep, in this case 1h.
-
+
ds_h = ds.interp_time(3600)
ds_h
@@ -412,14 +412,14 @@ Interpola
And to store the interpolated data in a new file.
-
+
ds_h.to_dfs("waves_3h.dfs2")
Interpolate to time axis of another dataset
Read some non-equidistant data typically found in observed data.
-
+
ts = mikeio.read("../data/waves.dfs0")
ts
@@ -434,10 +434,10 @@
+
dsi = ds.interp_time(ts)
-
+
dsi.time
DatetimeIndex(['2004-01-01 01:00:00', '2004-01-01 02:00:00',
@@ -455,13 +455,13 @@
+
dsi["Sign. Wave Height"].shape
(24, 31, 31)
-
+
ax = dsi["Sign. Wave Height"].sel(x=250, y=1200).plot(marker='+')
ts["Sign. Wave Height"].plot(ax=ax,marker='+')
@@ -479,7 +479,7 @@ Model validation
In the example below we calculate this metric using the model data interpolated to the observed times.
For a more elaborate model validation library which takes care of these things for you as well as calculating a number of relevant metrics, take a look at `ModelSkill.
Use np.nanmean to skip NaN.
-
+
ts["Sign. Wave Height"]
<mikeio.DataArray>
@@ -490,7 +490,7 @@ Model validation
values: [0.06521, 0.06771, ..., 0.0576]
-
+
dsi["Sign. Wave Height"].sel(x=250, y=1200)
<mikeio.DataArray>
@@ -501,7 +501,7 @@ Model validation
values: [0.0387, 0.03939, ..., nan]
-
+
diff = (ts["Sign. Wave Height"] - dsi["Sign. Wave Height"].sel(x=250, y=1200))
diff.plot()
@@ -512,7 +512,7 @@ Model validation
-
+
mae = np.abs(diff).nanmean().to_numpy()
mae
@@ -522,7 +522,7 @@ Model validation
Clean up
-
+
import os
os.remove("waves_3h.dfs2")
diff --git a/examples/dfs2/bathy.html b/examples/dfs2/bathy.html
index 5960c44c9..a1f61fc0f 100644
--- a/examples/dfs2/bathy.html
+++ b/examples/dfs2/bathy.html
@@ -266,7 +266,7 @@ Dfs2 - Bathymetric data
@@ -377,11 +377,11 @@ Dfs2 - Bathymetric data
GEBCO Compilation Group (2020) GEBCO 2020 Grid (doi:10.5285/a29c5465-b138-234d-e053-6c86abc040b9)
-
+
import xarray
import mikeio
-
+
ds = xarray.open_dataset("../../data/gebco_2020_n56.3_s55.2_w12.2_e13.1.nc")
ds
@@ -753,7 +753,7 @@ Dfs2 - Bathymetric data
Coordinates: (2)
Data variables:
elevation (lat, lon) int16 114kB ...
-Attributes: (8)- Conventions :
- CF-1.6
- title :
- The GEBCO_2020 Grid - a continuous terrain model for oceans and land at 15 arc-second intervals
- institution :
- On behalf of the General Bathymetric Chart of the Oceans (GEBCO), the data are held at the British Oceanographic Data Centre (BODC).
- source :
- The GEBCO_2020 Grid is the latest global bathymetric product released by the General Bathymetric Chart of the Oceans (GEBCO) and has been developed through the Nippon Foundation-GEBCO Seabed 2030 Project. This is a collaborative project between the Nippon Foundation of Japan and GEBCO. The Seabed 2030 Project aims to bring together all available bathymetric data to produce the definitive map of the world ocean floor and make it available to all.
- history :
- Information on the development of the data set and the source data sets included in the grid can be found in the data set documentation available from https://www.gebco.net
- references :
- DOI: 10.5285/a29c5465-b138-234d-e053-6c86abc040b9
- comment :
- The data in the GEBCO_2020 Grid should not be used for navigation or any purpose relating to safety at sea.
- node_offset :
- 1.0
ds.elevation.plot();
@@ -784,7 +784,7 @@ Dfs2 - Bathymetric data
-
+
ds.elevation.sel(lon=12.74792, lat=55.865, method="nearest")
+Attributes: (7)
Check ordering of dimensions, should be (y,x)
-
+
ds.elevation.dims
('lat', 'lon')
-
+
el = ds.elevation.values
el.shape
@@ -1171,37 +1171,37 @@ Dfs2 - Bathymetric data
Check that axes are increasing, S->N W->E
-
+
ds.lat.values[0],ds.lat.values[-1]
(55.20208333333332, 56.29791666666665)
-
+
ds.lat.values[0] < ds.lat.values[-1]
True
-
+
ds.lon.values[0],ds.lon.values[-1]
(12.20208333333332, 13.097916666666663)
-
+
el[0,0] # Bottom left
-8
-
+
el[-1,0] # Top Left
-31
-
+
geometry = mikeio.Grid2D(x=ds.lon.values, y=ds.lat.values, projection="LONG/LAT")
geometry
@@ -1211,7 +1211,7 @@ Dfs2 - Bathymetric data
projection: LONG/LAT
-
+
da = mikeio.DataArray(data=el,
item=mikeio.ItemInfo("Elevation", mikeio.EUMType.Total_Water_Depth),
geometry=geometry,
@@ -1226,7 +1226,7 @@ Dfs2 - Bathymetric data
geometry: Grid2D (ny=264, nx=216)
-
+
da.plot();
@@ -1236,7 +1236,7 @@ Dfs2 - Bathymetric data
-
+
da.plot(cmap='coolwarm', vmin=-100, vmax=100);
@@ -1246,10 +1246,10 @@ Dfs2 - Bathymetric data
-
+
da.to_dfs("gebco.dfs2")
-
+
ds = mikeio.read("gebco.dfs2")
ds.Elevation.plot()
@@ -1262,7 +1262,7 @@ Dfs2 - Bathymetric data
Clean up
-
+
import os
os.remove("gebco.dfs2")
diff --git a/examples/dfs2/gfs.html b/examples/dfs2/gfs.html
index 3c0467a6e..461e0bad8 100644
--- a/examples/dfs2/gfs.html
+++ b/examples/dfs2/gfs.html
@@ -266,7 +266,7 @@ Dfs2 - Meteo data
@@ -380,13 +380,13 @@ Dfs2 - Meteo data
-
+
import xarray
import pandas as pd
import mikeio
The file gfs_wind.nc contains a small sample of the GFS forecast data downloaded via their OpenDAP service
-
+
ds = xarray.open_dataset('../../data/gfs_wind.nc')
ds
@@ -760,30 +760,30 @@ Dfs2 - Meteo data
msletmsl (time, lat, lon) float32 10kB ...
ugrd10m (time, lat, lon) float32 10kB ...
vgrd10m (time, lat, lon) float32 10kB ...
-Attributes: (4)
Running a Mike 21 HD model, needs at least three variables of meteorological forcing * Mean Sea Level Pressure * U 10m * V 10m
Let’s take a look the U 10m
-
+
ds.ugrd10m.isel(time=0).plot();
@@ -797,7 +797,7 @@ Dfs2 - Meteo data
Convert to dfs2
Time
-
+
time = pd.DatetimeIndex(ds.time)
time
@@ -809,36 +809,36 @@ Time
Variable types
-
+
mikeio.EUMType.Air_Pressure
Air Pressure
-
+
mikeio.EUMType.Air_Pressure.units
[hectopascal, millibar]
-
+
mikeio.EUMType.Wind_Velocity
Wind Velocity
-
+
mikeio.EUMType.Wind_Velocity.units
[meter per sec, feet per sec, miles per hour, km per hour, knot]
-
+
mslp = ds.msletmsl.values / 100 # conversion from Pa to hPa
u = ds.ugrd10m.values
v = ds.vgrd10m.values
-
+
geometry = mikeio.Grid2D(x=ds.lon.values, y=ds.lat.values, projection="LONG/LAT")
geometry
@@ -848,14 +848,14 @@ Variable types
projection: LONG/LAT
-
+
from mikeio import ItemInfo, EUMType, EUMUnit
mslp_da = mikeio.DataArray(data=mslp,time=time, geometry=geometry, item=ItemInfo("Mean Sea Level Pressure", EUMType.Air_Pressure, EUMUnit.hectopascal))
u_da = mikeio.DataArray(data=u,time=time, geometry=geometry, item=ItemInfo("Wind U", EUMType.Wind_Velocity, EUMUnit.meter_per_sec))
v_da = mikeio.DataArray(data=v,time=time, geometry=geometry, item=ItemInfo("Wind V", EUMType.Wind_Velocity, EUMUnit.meter_per_sec))
-
+
mds = mikeio.Dataset([mslp_da, u_da, v_da])
mds
@@ -869,11 +869,11 @@ Variable types
2: Wind V <Wind Velocity> (meter per sec)
-
+
mds.to_dfs("gfs.dfs2")
Clean up
-
+
import os
os.remove("gfs.dfs2")
diff --git a/examples/dfs2/index.html b/examples/dfs2/index.html
index 580c5acd0..c5ff0ca9d 100644
--- a/examples/dfs2/index.html
+++ b/examples/dfs2/index.html
@@ -228,7 +228,7 @@ Dfs2 examples
diff --git a/examples/index.html b/examples/index.html
index b21302007..3cc800c42 100644
--- a/examples/index.html
+++ b/examples/index.html
@@ -288,7 +288,7 @@ Examples
@@ -420,7 +420,7 @@ Examples
-
+
Dfs2 - Bathymetric data
@@ -428,7 +428,7 @@ Examples
Convert GEBCO 2020 NetCDF to dfs2
-
+
Dfs2 - Meteo data
@@ -436,7 +436,7 @@ Examples
Conversion of NetCDF from Global Forecasting System to Dfs2
-
+
Dfs2 examples
@@ -447,7 +447,7 @@ Examples
-
+
Dfsu - 2D interpolation
@@ -455,7 +455,7 @@ Examples
Interpolate dfsu data to a grid, save as dfs2 and geotiff. Interpolate dfsu data to another mesh.
-
+
Generic dfs processing
@@ -463,7 +463,7 @@ Examples
Tools and methods that applies to any type of dfs files.
-
+
Time interpolation
@@ -866,8 +866,8 @@ Examples