This notebook demonstrates how to interpolate dfsu data to a grid, how to save the gridded data as dfs2 and geotiff. It also shows how to interpolate dfsu data to another mesh.
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
from mikeio import Dfsu, Mesh, Dfs2
from mikeio.spatial import Grid2D
C:\Users\jem\AppData\Local\Temp/ipykernel_10736/593968638.py:4: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()` set_matplotlib_formats('png')
dfs = Dfsu('../tests/testdata/wind_north_sea.dfsu')
ds = dfs.read(items=['Wind speed'])
dfs.plot(figsize=(8,6));
A Grid2d can be created either by giving x- and y-vectors, or by giving a bounding box and shape or spacing. The Grid2D can saved to mesh.
x = np.linspace(0,10,6)
y = [0,1.0,2.0]
g = Grid2D(x,y)
g
<mikeio.Grid2D> x-axis: nx=6 points from x0=0 to x1=10 with dx=2 y-axis: ny=3 points from y0=0 to y1=2 with dy=1 Number of grid points: 18
g.to_mesh('gridmesh.mesh')
from mikeio import Mesh
Mesh('gridmesh.mesh').plot(plot_type='mesh_only')
<AxesSubplot:>
g = dfs.get_overset_grid(shape=(50,40))
g
<mikeio.Grid2D> x-axis: nx=50 points from x0=-1.47424 to x1=8.74809 with dx=0.208619 y-axis: ny=40 points from y0=49.9408 to y1=55.259 with dy=0.136363 Number of grid points: 2000
g.bbox
BoundingBox(left=-1.5785510642197123, bottom=49.87266744109285, right=8.85240247077194, top=55.327192557197655)
g.x0
-1.4742415288697956
interpolant = dfs.get_2d_interpolant(g.xy, n_nearest=1)
dsi = dfs.interp2d(ds, *interpolant, shape=(g.ny, g.nx))
dsi
<mikeio.Dataset> Dimensions: (6, 40, 50) Time: 2017-10-27 00:00:00 - 2017-10-27 05:00:00 Items: 0: Wind speed <Wind speed> (meter per sec)
datgrid = dsi["Wind speed"][0]
plt.pcolormesh(g.xx, g.yy, datgrid);
C:\Users\jem\AppData\Local\Temp/ipykernel_10736/1104355596.py:2: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later. plt.pcolormesh(g.xx, g.yy, datgrid);
dat_orig = ds[0][0,:]
dfs.plot(dat_orig, show_mesh=False, figsize=(9,7));
import rasterio
from rasterio.transform import from_origin
# Dcoumentation https://rasterio.readthedocs.io/en/latest/index.html
with rasterio.open(
'wind.tif',
'w',
driver='GTiff',
height=g.ny,
width=g.nx,
count=1,
dtype=datgrid.dtype,
crs='+proj=latlong',
transform=from_origin(g.bbox.left, g.bbox.top, g.dx, g.dy)
) as dst:
dst.write(np.flipud(datgrid), 1)
Note: data needs to be flipped upside down
dsi.flipud()
<mikeio.Dataset> Dimensions: (6, 40, 50) Time: 2017-10-27 00:00:00 - 2017-10-27 05:00:00 Items: 0: Wind speed <Wind speed> (meter per sec)
dfs2_filename = 'wind_north_sea_interpolated.dfs2'
dfs2 = Dfs2()
coordinate = [dfs.projection_string, g.x0, g.y0, 0]
dfs2.write(dfs2_filename, data=dsi, coordinate=coordinate, dx=g.dx, dy=g.dy)
Interpolate the data from this coarse mesh onto a finer resolution mesh
msh = Mesh('../tests/testdata/north_sea_2.mesh')
dfs.n_elements, msh.n_elements # new mesh has more elements
(958, 2259)
xy = msh.element_coordinates[:,:2]
interpolant = dfs.get_2d_interpolant(xy, n_nearest=4)
dsi = dfs.interp2d(ds, *interpolant)
dat_orig = ds[0][0,:]
dfs.plot(dat_orig, figsize=(9,7));
dat_interp = dsi[0][0,:]
msh.plot(dat_interp, figsize=(9,7));
nan_elements = msh.element_ids[np.isnan(dat_interp)]
nan_elements
array([ 249, 451, 1546])
dfs.contains(msh.element_coordinates[nan_elements,:2])
array([False, False, False])
interpolant = dfs.get_2d_interpolant(xy, n_nearest=4, extrapolate=True)
dat_interp = dfs.interp2d(dat_orig, *interpolant)
n_nan_elements = sum(np.isnan(dat_interp))
n_nan_elements
0
We want to interpolate scatter data onto an existing mesh and create a new dfsu with the interpolated data.
from mikeio.spatial import dist_in_meters
from mikeio.interpolation import get_idw_interpolant
dfs = Dfsu('../tests/testdata/wind_north_sea.dfsu')
dfs.plot(plot_type="mesh_only");
# 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
array([[ 1. , 50. , 1. ], [ 4. , 52. , 3. ], [ 8. , 55. , 2. ], [-1. , 55. , 1.5]])
Let's first try the approx for a single element:
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)
dfs.plot(dati, title="Interpolated scatter data");
dfs.write("interpolated_scatter.dfsu", [dati], start_time=dfs.start_time)
import os
os.remove("wind_north_sea_interpolated.dfs2")
os.remove("gridmesh.mesh")
os.remove("wind.tif")
os.remove("interpolated_scatter.dfsu")