#!/usr/bin/env python
# coding: utf-8

# # Global Forecasting System - Meteorological forecast

# In[1]:


import xarray
import pandas as pd
import mikeio


# The file `gfs_wind.nc` contains a small sample of the [GFS](https://nomads.ncep.noaa.gov/) forecast data downloaded via their OpenDAP service

# In[2]:


ds = xarray.open_dataset('../tests/testdata/gfs_wind.nc')


# Running a Mike 21 HD model, needs at least three variables of meteorological forcing
# * Mean Sea Level Pressure
# * U 10m
# * V 10m

# In[3]:


ds.msletmsl


# In[4]:


ds.ugrd10m


# In[5]:


ds.vgrd10m


# In[6]:


ds.ugrd10m[0].plot();


# ## Convert to dfs2

# ## Time

# In[7]:


time = pd.DatetimeIndex(ds.time)
time


# ## Variable types

# In[8]:


mikeio.EUMType.Air_Pressure


# In[9]:


mikeio.EUMType.Air_Pressure.units


# In[10]:


mikeio.EUMType.Wind_Velocity


# In[11]:


mikeio.EUMType.Wind_Velocity.units


# In[12]:


mslp = ds.msletmsl.values / 100 # conversion from Pa to hPa
u = ds.ugrd10m.values
v = ds.vgrd10m.values


# In[13]:


geometry = mikeio.Grid2D(x=ds.lon.values, y=ds.lat.values, projection="LONG/LAT")
geometry


# In[14]:


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))


# In[15]:


mds = mikeio.Dataset([mslp_da, u_da, v_da])
mds


# In[16]:


mds.to_dfs("gfs.dfs2")


# ## Clean up (don't run this)

# In[17]:


import os
os.remove("gfs.dfs2")