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 mikeio
ds = mikeio.read("../tests/testdata/wind_north_sea.dfsu", items="Wind speed")
ds
<mikeio.Dataset> dims: (time:6, element:958) time: 2017-10-27 00:00:00 - 2017-10-27 05:00:00 (6 records) geometry: Dfsu2D (958 elements, 570 nodes) items: 0: Wind speed <Wind speed> (meter per sec)
da = ds.Wind_speed
da.plot();
g = da.geometry.get_overset_grid(dx=0.1)
g
<mikeio.Grid2D> x: [-1.563, -1.463, ..., 8.837] (nx=105, dx=0.1) y: [49.9, 50, ..., 55.3] (ny=55, dy=0.1) projection: LONG/LAT
da_grid = da.interp_like(g)
da_grid
<mikeio.DataArray> name: Wind speed dims: (time:6, y:55, x:105) time: 2017-10-27 00:00:00 - 2017-10-27 05:00:00 (6 records) geometry: Grid2D (ny=55, nx=105)
da_grid.plot();
da_grid.to_dfs("wind_north_sea_interpolated.dfs2")
xr_da = da_grid.to_xarray()
xr_da.to_netcdf("wind_north_sea_interpolated.nc")
Install rasterio by running this in a command prompt before running this notebook
$ conda install -c conda-forge rasterio
Or if you prefer to avoid conda, here is how: https://rasterio.readthedocs.io/en/latest/installation.html#windows
import numpy as np
import rasterio
from rasterio.transform import from_origin
# Dcoumentation https://rasterio.readthedocs.io/en/latest/index.html
with rasterio.open(
fp='wind.tif',
mode='w',
driver='GTiff',
height=g.ny,
width=g.nx,
count=1,
dtype=da.dtype,
crs='+proj=latlong', # adjust accordingly for projected coordinate systems
transform=from_origin(g.bbox.left, g.bbox.top, g.dx, g.dy)
) as dst:
dst.write(np.flipud(da_grid[0].to_numpy()), 1) # first time_step
Interpolate the data from this coarse mesh onto a finer resolution mesh
msh = mikeio.Mesh('../tests/testdata/north_sea_2.mesh')
msh
<Mesh> number of elements: 2259 number of nodes: 1296 projection: LONG/LAT
dsi = da.interp_like(msh)
dsi
<mikeio.DataArray> name: Wind speed dims: (time:6, element:2259) time: 2017-10-27 00:00:00 - 2017-10-27 05:00:00 (6 records) geometry: Dfsu2D (2259 elements, 1296 nodes)
da[0].plot(figsize=(9,7), show_mesh=True);
dsi[0].plot(figsize=(9,7), show_mesh=True);
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])
dat_interp = da.interp_like(msh, extrapolate=True)
n_nan_elements = np.sum(np.isnan(dat_interp.values))
n_nan_elements
0
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('../tests/testdata/wind_north_sea.dfsu')
dfs.geometry.plot.mesh();
# 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.geometry.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.geometry.n_elements))
for j in range(dfs.geometry.n_elements):
dist = dist_in_meters(scatter[:,:2], dfs.geometry.element_coordinates[j,:2])
w = get_idw_interpolant(dist, p=2)
dati[0,j] = np.dot(scatter[:,2], w)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[25], line 1 ----> 1 dati = np.zeros((1,dfs.n_elements)) 2 for j in range(dfs.n_elements): 3 dist = dist_in_meters(scatter[:,:2], dfs.geometry.element_coordinates[j,:2]) AttributeError: 'Dfsu2DH' object has no attribute 'n_elements'
da = mikeio.DataArray(data=dati, geometry=dfs.geometry, time=dfs.start_time)
da
da.plot(title="Interpolated scatter data");
da.to_dfs("interpolated_scatter.dfsu")
import os
os.remove("wind_north_sea_interpolated.dfs2")
os.remove("wind_north_sea_interpolated.nc")
os.remove("wind.tif")
os.remove("interpolated_scatter.dfsu")