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

# # Dfsu - Export to shapefile
# 
# 1. Read selected item and timestep from dfsu
# 2. Extract geometry
# 3. Create GeoPandas dataframe
# 4. Save to ESRI shapefile

# In[1]:


# conda install geopandas


# ## Step 1. read the selected data

# In[2]:


from mikeio import Dfsu

dfsu = Dfsu("../tests/testdata/wind_north_sea.dfsu")
ds = dfsu.read()
ws = ds["Wind speed"][0]
dfsu.plot(ws, label="Wind speed")


# ## Step 2. extract geometry

# In[3]:


shp = dfsu.to_shapely()
type(shp)


# Geopandas does not like multipolygon, it should be a list of polygons

# In[4]:


poly_list = [p for p in shp]


# ## Step 3. Create a geopandas dataframe

# In[5]:


import pandas as pd
import geopandas as gpd

df = pd.DataFrame({'wind_speed':ws})
df.head()


# In[6]:


gdf = gpd.GeoDataFrame(df,geometry=poly_list)


# ## Step 4. Save to shapefile

# In[7]:


gdf.to_file("wind_speed.shp")


# ## Step 5...
# Do further work in QGIS
# 
# ![QGIS](../images/dfsu_qgis.png)

# Would you prefer to have this workflow to be a method on the `mikeio.Dfsu` class?
# 
# Post an issue on [GitHub](https://github.com/DHI/mikeio/issues) !

# # Contour lines

# In[8]:


import matplotlib.pyplot as plt


# In[9]:


# get coordinates
ec = dfsu.element_coordinates
lon = ec[:,0]
lat = ec[:,1]


# In[10]:


# Select item and timestep
m = ds['Wind speed'][0]


# In[11]:


# Interpolate to cartesian grid
import numpy as np
from  scipy.interpolate import griddata
numcols, numrows = 200, 200
xi = np.linspace(lon.min(), lon.max(), numcols)
yi = np.linspace(lat.min(), lat.max(), numrows)
xi, yi = np.meshgrid(xi, yi)

grid_z = griddata(points=ec[:,0:2],values=m,xi=(xi,yi),method='cubic')


# In[12]:


contour_levels=np.arange(4, 14, 0.5)
cn = plt.contour(xi,yi,grid_z,levels=contour_levels)


# In[13]:


from shapely.geometry import LineString

poly_list = []

for i in range(len(cn.collections)):
    p = cn.collections[i].get_paths()[0]
    v = p.vertices
    x = v[:,0]
    y = v[:,1]
    poly = LineString([(i[0], i[1]) for i in zip(x,y)])
    if(poly.is_empty):
        print(f"{i} is empty")
    poly_list.append(poly)


# In[14]:


# Clip to domain
domain = dfsu.to_shapely().buffer(0)
poly_list = [p.intersection(domain) for p in poly_list]


# In[15]:


# Create GeoDataframe
df = pd.DataFrame({'wind_speed':contour_levels})
gdf = gpd.GeoDataFrame(df,geometry=poly_list)
gdf.head()


# In[16]:


# export shapefile
gdf.to_file("wind_speed_contours.shp")


# ![QGIS](../images/dfsu_qgis_contours.png)

# # Clean up

# In[17]:


import os

files = ["wind_speed","wind_speed_contours"]

exts = ["cpg","dbf","shp","shx"]

for file in files:
    for ext in exts:
        filename = f"{file}.{ext}"
        if os.path.exists(filename): os.remove(filename)


#