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

# # Mesh
# * read mesh file
# * plot mesh 
# * convert to shapely
# * check if point is inside or outside mesh
# * subset mesh, plot subset
# * change z values
# * change boundary codes

# In[35]:


import matplotlib.pyplot as plt
from matplotlib_inline.backend_inline import set_matplotlib_formats
set_matplotlib_formats('png')
plt.rcParams["figure.figsize"] = (6,6)

import mikeio


# In[36]:


msh = mikeio.Mesh("../tests/testdata/odense_rough.mesh")
msh


# In[37]:


msh.plot()
msh.plot.boundary_nodes(boundary_names=['Land','Open boundary']);


# # Convert mesh to shapely
# Convert mesh to [shapely](https://shapely.readthedocs.io/en/latest/manual.html) MultiPolygon object, requires that the `shapely` library is installed.

# In[38]:


mp = msh.to_shapely()
mp


# Now a lot of methods are available

# In[39]:


mp.area


# In[40]:


mp.bounds


# In[41]:


domain = mp.buffer(0)
domain


# In[42]:


open_water = domain.buffer(-500)

coastalzone = domain - open_water
coastalzone


# Find if points are inside the domain

# In[43]:


from shapely.geometry import Point

p1 = Point(216000, 6162000)
p2 = Point(220000, 6156000)
print(mp.contains(p1))
print(mp.contains(p2))


# ## Mesh class can also check if a mesh contains points 

# In[44]:


p1p2 = [[216000, 6162000], [220000, 6156000]]
msh.geometry.contains(p1p2)


# In[45]:


ax = msh.plot()
ax.scatter(p1.x, p1.y, marker="*", s=200, c="red", label="inside")
ax.scatter(p2.x, p2.y, marker="+", s=200, c="green", label="outside")
ax.legend();


# # Subset mesh
# Select only elements with more than 3m depth. Plot these elements. 

# In[46]:


zc = msh.element_coordinates[:,2]


# # Change z values and boundary code
# Assume that we want to have a minimum depth of 2 meters and change the open boundary (code 2) to a closed one (code 1). 

# In[47]:


print(f'max z before: {msh.node_coordinates[:,2].max()}')
zc = msh.node_coordinates[:,2]
zc[zc>-2] = -2
msh.zn = zc
print(f'max z after: {msh.node_coordinates[:,2].max()}')


# In[48]:


c = msh.geometry.codes
c[c==2] = 1
msh.geometry.codes = c


# In[49]:


msh.geometry.codes