<jupyter_start><jupyter_text>Input data Weather forcingThe first example uses data from ECMWF ERA5 obtained from [CDS](https://cds.climate.copernicus.eu/!/home)A common format for multidimensional data is [NetCDF](https://www.unidata.ucar.edu/software/netcdf/).An excellent python package for working with NetCDF in Python is [Xarray](http://xarray.pydata.org/en/stable/). To work with NetCDF data first we need to install the `netcdf4` and `xarray` packages<jupyter_code>! pip install netcdf4 xarray<jupyter_output><empty_output><jupyter_text>Opening a NetCDF file with Xarray returns an `xarray.Dataset`<jupyter_code>import xarray as xr

ds = xr.open_dataset("data/era5_small.nc")
ds<jupyter_output><empty_output><jupyter_text>One of the datam variables is the t2m (Air temperature at 2m above the ground) which we can access in two different ways:1. Using the string "t2m"2. Or as a property of the Dataset object `.t2m`<jupyter_code>ds["t2m"]
ds.t2m<jupyter_output><empty_output><jupyter_text>To select the first timestep, we can use the `xarray.Dataset.isel` method.<jupyter_code>ds0 = ds.isel(time=0)
ds0<jupyter_output><empty_output><jupyter_text>The dataset from the first timestep is also a Dataset, except it no longer has a time dimension.To plot a map of the air temperature we can use select the `t2m` property and call the plot method.<jupyter_code>ds0.t2m.plot()<jupyter_output><empty_output><jupyter_text>Pressure Atmospheric pressure comes in different flavors.`Surface pressure` (sp) is the pressure at ground level, i.e. following the terrain.In the following plot you can clearly see the decrease in surface pressure with altitude in areas with mountains.This is *NOT* the variable used by MIKE 21 FM HD.<jupyter_code>ds0.sp.plot()<jupyter_output><empty_output><jupyter_text>MIKE 21 FM HD expects the atmospheric pressure reduced to Mean Sea Level. (`msl`) In the following plot the pressure field has a much smaller range, since the variation with altitude have been removed.<jupyter_code>ds0.msl.plot()<jupyter_output><empty_output><jupyter_text>Wind The wind velocity components U,V at 10m above the ground, are found in the two data variables, `u10` and `v10`<jupyter_code>ds0.u10.long_name
ds0.v10.long_name<jupyter_output><empty_output><jupyter_text>Convert to dfs2  Sub-regionIn this case the NetCDF file covers a much larger area, than is required for our modelling work, so first step is to create a spatial subset to match our Area of Interest.Since the latitude dimension is ordered in decreasing order in this dataset, the slice has to go from high to low.<jupyter_code>ds_aoi = ds.sel(longitude=slice(10,30), 
                latitude=slice(60,45)) # N -> S
ds_aoi<jupyter_output><empty_output><jupyter_text>Make a plot to check that the area is the one we expect.<jupyter_code>ds_aoi.isel(time=0).u10.plot(cmap='jet')<jupyter_output><empty_output><jupyter_text>TimeNow we are ready to convert this data to dfs2.First step is to create a time axis, which can be understood by MIKE IO.<jupyter_code>import pandas as pd

time = pd.DatetimeIndex(ds.time)
time<jupyter_output><empty_output><jupyter_text>The second step is to collect the variables that we need, and at the same time do unit conversion where it is necessary. In this case for pressure and temperature.<jupyter_code>data = [ds_aoi.msl.values / 100.0, # Pa -> hPa
        ds_aoi.u10.values,
        ds_aoi.v10.values,
        ds_aoi.t2m.values - 273.15] # K -> °C<jupyter_output><empty_output><jupyter_text>The third step is to create the proper EUM item information for each variable.<jupyter_code>from mikeio.eum import EUMType, EUMUnit, ItemInfo

items = [
        ItemInfo("MSLP", EUMType.Pressure, EUMUnit.hectopascal),
        ItemInfo("U 10m", EUMType.Wind_Velocity, EUMUnit.meter_per_sec),
        ItemInfo("V 10m", EUMType.Wind_Velocity, EUMUnit.meter_per_sec),
        ItemInfo("Temperature", EUMType.Temperature, EUMUnit.degree_Celsius)
]<jupyter_output><empty_output><jupyter_text>The fourth step is to gather the data, time and item information in a `mikeio.Dataset`<jupyter_code>from mikeio import Dataset

my_ds = Dataset(data=data, time=time, items=items)
my_ds
ds_aoi.longitude.values[0], ds.longitude.values[-1]
dx = ds_aoi.longitude.values[1] - ds_aoi.longitude.values[0]
dx
ds_aoi.latitude.values[0], ds_aoi.latitude.values[-1]
dy = ds_aoi.latitude.values[0] - ds_aoi.latitude.values[1] # dx, dy should be positive!
dy
type(data[0])<jupyter_output><empty_output><jupyter_text>The dimensions for at  DFS2 are expected to be (t, y, x), which matches the ones used by this dataset.This convention is recommended by the [CF convention](http://cfconventions.org/cf-conventions/cf-conventions.htmldimensions)<jupyter_code>ds.Conventions
ds.u10.dims
data[0].shape
lon_origin = ds_aoi.longitude.values[0]
lon_origin
lat_origin = ds_aoi.latitude.values[-1] # Southern-most point
lat_origin
# This dataset North-South orientation matches the expected, but for other dataset you might need to flip the data upside-down
# It can be done like this
# my_ds = my_ds.flipud()<jupyter_output><empty_output><jupyter_text>The final step is to write the dataset to a dfs2 file, where we also provide information on the coordinate reference system to use (i.e. projection, origin, spacing).<jupyter_code>from mikeio import Dfs2

orientation = 0.0

dfs = Dfs2()
dfs.write("era5_aoi.dfs2", data=my_ds, coordinate=["LONG/LAT", lat_origin, lon_origin, orientation], dx=dx, dy=dy)<jupyter_output><empty_output><jupyter_text>Screenshot of U 10m from MIKE Zero![](images/era5_u10_dfs2.png)  Unstructured grid forcingE.g. roughness map<jupyter_code>from mikeio import Mesh, Dfsu

msh = Mesh("data/FakeLake.dfsu")
msh.plot()
import numpy as np
import matplotlib.pyplot as plt

z = np.linspace(0,-30)
manning = 32 + (40 - 32)*(1 - np.exp(z/10))

plt.plot(z, manning)
plt.xlim(0,-30)
plt.ylabel("Manning")
plt.xlabel("Depth (m)")
plt.title("This is only an example...")<jupyter_output><empty_output><jupyter_text>Extract the Z values in each element.<jupyter_code>ze = msh.element_coordinates[:,2]<jupyter_output><empty_output><jupyter_text>Create a new array of manning numbers based on the depth. This is only a hypothetical example of a data transformation based on the data in the file, not a best practise for bed friction.<jupyter_code>me = 32 + (40 - 32)*(1 - np.exp(ze/10))<jupyter_output><empty_output><jupyter_text>The primary use case of DFS files is for time-varying data, so in order to use this time-independent data, we have to introduce a new axis, a time dimension.One simple way to do this is to use the numpy function `atleast_2d`.<jupyter_code>me2d = np.atleast_2d(me) # Introduce time axis
me2d.shape<jupyter_output><empty_output><jupyter_text>Now the array is a 2d array<jupyter_code>from mikeio.eum import ItemInfo
from mikeio import Dataset

my_ds = Dataset(data=[me2d],
                time="2000-1-1", # For static data, this is not used by MIKE 21
                items=[ItemInfo("Roughness", EUMType.Mannings_M)])
my_ds
dfs = Dfsu("data/FakeLake.dfsu") # use this file as template
dfs.write("manning.dfsu", data=my_ds)
dfs = Dfsu("manning.dfsu")
dfs
ds = dfs.read()
dfs.plot(ds["Roughness"][0], title="Manning")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1, 'fraction_numerical': 0.018620689655172412, 'mean_word_length': 5.5501355013550135, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 88, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '514463', 'license': 'MIT', 'n_tokens_mistral': 2609, 'n_tokens_neox': 2497, 'n_words': 747, 'path': 'mini_book/input_data.ipynb', 'repo_name': 'DHI/book-learn-mikeio-fmskill', 'size': '7250'}