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

# # Show FESOM mesh

# This notebook shows how to visualize FESOM grid. We are going to use *%matplotlib notebook* option in Jupyter notebooks in order to make the plots interactive.

# In[2]:


import sys
sys.path.append("../")

import pyfesom as pf
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
import numpy as np
get_ipython().run_line_magic('matplotlib', 'notebook')
from matplotlib import cm


# First we have to load the FESOM mesh. We only need information about surface nodes and elements (triangles), so we put *get3d* to False.

# In[3]:


meshpath  ='../../../../FESOM/mesh/'
mesh = pf.fesom_mesh(meshpath, get3d=False)


# You can have a look at the basic information about the mesh. Note that the information about 3d nods are not loaded.

# In[4]:


mesh


# Here we remove cyclic elements from the variables that contain information about elements and volumes (2d volumes, really areas).

# In[5]:


elem2=mesh.elem[mesh.no_cyclic_elem,:]
voltri = mesh.voltri[mesh.no_cyclic_elem]


# Setup a map and convert positions of the elements to the map projection

# In[19]:


map = Basemap(projection='robin',lon_0=0, resolution='l')
x, y = map(mesh.x2, mesh.y2)


# Plot triangular mesh:
# 
#     x     - lon (x) coordinates
#     y     - lat (y) coordinates
#     elem2 - 3d array, where each raw have trhee indexes of the x and y. So, each raw describes triangle 

# In[20]:


plt.figure(figsize=(9,5))
map.drawmapboundary(fill_color='0.9')
map.drawcoastlines()
plt.triplot(x, y, elem2, lw=0.2,color='k');


# The resulting plot is interactive. You can zoom to specific areas and do some other basic manupulations by pressing icons on the lower left. 
# 
# IMPORTANT - when the figure is in the interactive mode, it considered to be the active one, so the output of the following ploting functions will apear on this figure as well. To avoid this when you finish to work with the figure, you have to press the "switch off" button at the upper right of the figure.

# One can plot color codded volume of the triangles as well. This is is much more computationally expensive. 

# In[19]:


#this cell can be executed only once if you don't plan to change the Basemap settings 
map = Basemap(projection='robin',lon_0=0, resolution='l')
x, y = map(mesh.x2, mesh.y2)


# In[20]:


map.drawmapboundary(fill_color='0.9')
map.drawcoastlines()
plt.tripcolor(x, y, elem2, \
              facecolors=mesh.voltri[mesh.no_cyclic_elem], \
              edgecolors='k',\
             cmap=cm.Spectral_r)
plt.colorbar(orientation='horizontal', pad=0.03);


# The Basemap can use [many different projections](http://matplotlib.org/basemap/users/mapsetup.html) and one of course can define the region to be plotted. Please refer to the [Basemap documentation](http://matplotlib.org/basemap/index.html) for details.

# ##  North Pole stereographic projection

# In[16]:


map = Basemap(projection='npstere',boundinglat=60,lon_0=0,resolution='l')
x, y = map(mesh.x2, mesh.y2)


# In[18]:


map.drawmapboundary(fill_color='0.9')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
#note that here we use *mesh.elem* directly since the projection can adequatelly process it.
im=plt.triplot(x, y, mesh.elem, lw=0.2,color='k')


# ## Mercator projection with etopo topography

# In[29]:


map = Basemap(projection='merc',llcrnrlat=30,urcrnrlat=70,\
            llcrnrlon=-30,urcrnrlon=30,lat_ts=20,resolution='l')

x, y = map(mesh.x2, mesh.y2)


# In[30]:


plt.figure(figsize=(8,9))
map.drawmapboundary(fill_color='0.9')
map.drawcoastlines()
plt.triplot(x, y, elem2, lw=0.2,color='k');
map.etopo()


# In[ ]: