from mikeio import Dfsu
filename = "../tests/testdata/HD2D.dfsu"
dfs = Dfsu(filename)
dfs
Dfsu2D Number of elements: 884 Number of nodes: 529 Projection: UTM-29 Items: 0: Surface elevation <Surface Elevation> (meter) 1: U velocity <u velocity component> (meter per sec) 2: V velocity <v velocity component> (meter per sec) 3: Current speed <Current Speed> (meter per sec) Time: 9 steps with dt=9000.0s 1985-08-06 07:00:00 -- 1985-08-07 03:00:00
ds = dfs.read(["Surface elevation","Current speed"]) # to read some variables
ds
<mikeio.Dataset> Dimensions: (9, 884) Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 Items: 0: Surface elevation <Surface Elevation> (meter) 1: Current speed <Current Speed> (meter per sec)
Dfsu.show_progress = True # Turn on progress bar, useful to get progress for long-running tasks
ds = dfs.read()
100%|█████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 14745.60it/s]
ds.describe()
Surface elevation | U velocity | V velocity | Current speed | |
---|---|---|---|---|
count | 7956.000000 | 7956.000000 | 7956.000000 | 7956.000000 |
mean | 0.034152 | -0.027504 | 0.051639 | 0.073658 |
std | 0.534890 | 0.088108 | 0.092440 | 0.119606 |
min | -0.758378 | -1.280627 | -0.811804 | 0.000082 |
25% | -0.579662 | -0.022283 | 0.019642 | 0.029917 |
50% | 0.101636 | -0.011019 | 0.041748 | 0.046848 |
75% | 0.462013 | 0.004649 | 0.062974 | 0.070259 |
max | 1.172707 | 0.056265 | 1.518996 | 1.714988 |
Find which element is nearest to POI.
idx = dfs.find_nearest_elements(606200, 6905480)
Extract a subset of the dataset from this element. (Discrete values, no interpolation)
selds = ds.isel(idx=idx)
selds
<mikeio.Dataset> Dimensions: (9,) Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 Items: 0: Surface elevation <Surface Elevation> (meter) 1: U velocity <u velocity component> (meter per sec) 2: V velocity <v velocity component> (meter per sec) 3: Current speed <Current Speed> (meter per sec)
Convert to a dataframe, for convenience.
df = selds.to_dataframe()
df.head()
Surface elevation | U velocity | V velocity | Current speed | |
---|---|---|---|---|
1985-08-06 07:00:00 | 0.459460 | 0.006372 | -0.007143 | 0.009572 |
1985-08-06 09:30:00 | 0.806965 | 0.010517 | 0.003438 | 0.011064 |
1985-08-06 12:00:00 | 0.100285 | 0.011300 | 0.012926 | 0.017169 |
1985-08-06 14:30:00 | -0.727009 | 0.010402 | 0.010033 | 0.014452 |
1985-08-06 17:00:00 | -0.579541 | 0.007293 | -0.000902 | 0.007349 |
df.plot()
<AxesSubplot:>
Assume that we interested in these 3 points only
pt1 = (606200, 6905480)
pt2 = (606300, 6905410)
pt3 = (606400, 6905520)
pts_x = [pt1[0], pt2[0], pt3[0]]
pts_y = [pt1[1], pt2[1], pt3[1]]
elem_ids = dfs.find_nearest_elements(pts_x, pts_y)
We can use these element ids either when we select the data from the complete dataset using the method isel() as shown above or already when we read the data from file (particular useful for files larger than memory)
ds_pts = dfs.read(elements=elem_ids)
ds_pts
100%|█████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 13878.21it/s]
<mikeio.Dataset> Dimensions: (9, 3) Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 Items: 0: Surface elevation <Surface Elevation> (meter) 1: U velocity <u velocity component> (meter per sec) 2: V velocity <v velocity component> (meter per sec) 3: Current speed <Current Speed> (meter per sec)
Let's take the area North of y=6905480
yc = dfs.element_coordinates[:,1]
elem_ids = dfs.element_ids[yc>6905480]
And find the maximum average current speed in this area in the last time step
item_num = 1
print(ds.items[item_num])
subset = ds.data[1][:,elem_ids]
subset_timeavg = subset.mean(axis=0)
idx = subset_timeavg.argmax()
coords = dfs.element_coordinates[idx,0:2].round(1)
print(f'Max current speed in area is found in {coords} and is {subset_timeavg[idx]:.3f}m/s')
U velocity <u velocity component> (meter per sec) Max current speed in area is found in [ 606674.6 6905718.3] and is 0.039m/s
Let us save the time averaged subset to a dfsu file.
outfilename1 = "HD2D_north.dfsu"
data = []
data.append(subset_timeavg.reshape(1,-1))
items = ds.items[item_num]
dfs.write(outfilename1, data, items=[items], elements=elem_ids)
c:\users\jem\source\mikeio\mikeio\dfsu.py:2342: UserWarning: No start time supplied. Using start time from source: 1985-08-06 07:00:00 as start time. warnings.warn( 100%|██████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 8355.19it/s]
First inspect the source file:
ds = dfs.read()
ds
100%|██████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 7420.63it/s]
<mikeio.Dataset> Dimensions: (9, 884) Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00 Items: 0: Surface elevation <Surface Elevation> (meter) 1: U velocity <u velocity component> (meter per sec) 2: V velocity <v velocity component> (meter per sec) 3: Current speed <Current Speed> (meter per sec)
from mikeio.eum import ItemInfo, EUMType
from mikeio import Dataset
sourcefilename = filename
outfilename2 = "HD2D_selected.dfsu"
starttimestep = 4
time = ds.time[starttimestep:]
data = []
data.append(ds['U velocity'][starttimestep:,:])
data.append(ds['V velocity'][starttimestep:,:])
items = [ItemInfo("eastward_sea_water_velocity", EUMType.u_velocity_component),
ItemInfo("northward_sea_water_velocity",EUMType.v_velocity_component)]
newds = Dataset(data,time,items)
dfs.write(outfilename2, newds) # Note, this method was previosly named create
100%|██████████████████████████████████████████████████████████████████████████████████| 5/5 [00:00<00:00, 9991.20it/s]
Read the newly created file to verify the contents.
newdfs = Dfsu(outfilename2)
newds2 = newdfs.read()
newds2
100%|██████████████████████████████████████████████████████████████████████████████████| 5/5 [00:00<00:00, 9972.19it/s]
<mikeio.Dataset> Dimensions: (5, 884) Time: 1985-08-06 17:00:00 - 1985-08-07 03:00:00 Items: 0: eastward_sea_water_velocity <u velocity component> (meter per sec) 1: northward_sea_water_velocity <v velocity component> (meter per sec)
Don't you have the original mesh? No problem - you can re-create it from the dfsu file...
outmesh = 'mesh_from_HD2D.mesh'
dfs.to_mesh(outmesh)
import os
os.remove(outfilename1)
os.remove(outfilename2)
os.remove(outmesh)