#
# Before opening Python execute:
# module load python3/miniconda3
# source activate base_plus

# Importation of modules and library
import numpy as np                         # Good module for matrix and matrix operation
import matplotlib as mpl
import matplotlib.colors as colors
import os                                  # Used to convert png to other format
import copy
import netCDF4

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pylab as plt

########################## Lecture de la pression au niveau de la mer et conversion de Pa à hPa
chemin_donnees_msl = '/sca/compilers_and_tools/python/demos-python3/data/era5_msl_ll_2019032512_1h.nc4'
fichier_netcdf_msl = netCDF4.Dataset(chemin_donnees_msl)
champ_msl = (1.0/100.0)*fichier_netcdf_msl.variables['msl'][0,:,:] # 1 [hPa] = 1/100 * [Pa]
lat = fichier_netcdf_msl.variables['latitude'][:]
lon = fichier_netcdf_msl.variables['longitude'][:]

fichier_netcdf_msl.close()
#########################

########################## Lecture du vent en x
chemin_donnees_u10 = '/sca/compilers_and_tools/python/demos-python3/data/era5_u10_ll_2019032512_1h.nc4'
fichier_netcdf_u10 = netCDF4.Dataset(chemin_donnees_u10)
champ_u10 = fichier_netcdf_u10.variables['u10'][0,:,:] # en m/s
#lat_u10 = fichier_netcdf_u10.variables['latitude'][:]
#lon_u10 = fichier_netcdf_u10.variables['longitude'][:]

fichier_netcdf_u10.close()
#########################

########################## Lecture du vent en y
chemin_donnees_v10 = '/sca/compilers_and_tools/python/demos-python3/data/era5_v10_ll_2019032512_1h.nc4'
fichier_netcdf_v10 = netCDF4.Dataset(chemin_donnees_v10)
champ_v10 = fichier_netcdf_v10.variables['v10'][0,:,:] # en m/s
#lat_v10 = fichier_netcdf_v10.variables['latitude'][:]
#lon_v10 = fichier_netcdf_v10.variables['longitude'][:]

fichier_netcdf_v10.close()
#########################

## Titre du graphique
titre = "ERA5 25 mars 2019 12Z PNM et vents à 10 m" 

## Nom final du fichier de l'image et qualité
image_output_file='./demo_msl_vent.gif'
image_output_dpi=150

###################### Faire la carte
## Création de la figure et choix de la projection
fig = plt.figure(figsize=(8,6))
ax = plt.axes(projection=ccrs.Miller())

## Conversion de lat-lon de 1D à 2D
lon2d, lat2d = np.meshgrid(lon, lat)

## Limite de la carte (ici en coordonnées géographiques)
# Lorsque le crs n'est pas spécifié, on suppose que c'est en crs.Geodetic()
ax.set_extent([-47,27,51,81])

## Features géographiques
# ax.coastlines(resolution='50m');
ax.add_feature(cfeature.OCEAN.with_scale('50m'))      # couche ocean
ax.add_feature(cfeature.LAND.with_scale('50m'))       # couche land
#ax.add_feature(cfeature.LAKES.with_scale('50m'))      # couche lac (polygone)
ax.add_feature(cfeature.LAKES.with_scale('50m'),facecolor='none',edgecolor='black',linewidth=0.5)      # couche lac (bordure)
ax.add_feature(cfeature.BORDERS.with_scale('50m'),linewidth=0.5)    # couche frontieres
#ax.add_feature(cfeature.RIVERS.with_scale('50m'))     # couche rivières 
coast = cfeature.NaturalEarthFeature(category='physical', scale='50m', facecolor='none', name='coastline',linewidth=0.5) # ajout de la couche cotière 
ax.add_feature(coast, edgecolor='black')

#states_provinces = cfeature.NaturalEarthFeature(category='cultural',name='admin_1_states_provinces_lines',scale='10m',facecolor='none') # Provinces/états
#ax.add_feature(states_provinces, edgecolor='gray')

## Valeurs minimales et maximales du champ de pression
val_min = 1000.0
val_max = 1030.0
val_int = 5.0
data_levels = np.arange(val_min,val_max+val_int,val_int)

## Dessin du champ de pression
data_contour = ax.contour(lon2d,lat2d,champ_msl,vmin=val_min,vmax=val_max,transform=ccrs.PlateCarree(),levels=data_levels,colors='k', extend='max' )
plt.clabel(data_contour,  data_levels, inline=True, fmt='%1i', fontsize=10)

## Dessin des vecteurs de vents
# Pour garder seulement un vecteur sur 10
skip=(slice(None,None,10),slice(None,None,10))

# Dessin des vecteurs
data_quiver = ax.quiver(lon2d[skip],lat2d[skip],champ_u10[skip],champ_v10[skip],color='red', scale = 250,transform=ccrs.PlateCarree())

# Pour avoir un vecteur exemple de vitesse
qk = ax.quiverkey(data_quiver, 0.8, 0.96, 10, "10 m/s", labelpos='E', coordinates='figure')

## Titre du graphique
plt.title(titre, fontsize=10)

## Pour s'assurer que tout est dessiné
fig.canvas.draw() 

## Pour élliminer le plus possible de bordure blanche
plt.tight_layout()

## Sauvegarde de la figure
fig.savefig('python.png', dpi=image_output_dpi, format='png',bbox_inches='tight') # Most backends support png, pdf, ps, eps and svg
os.system ('convert python.png ' + image_output_file)         # Convert python.png to 'image_output_file', i.e. .gif
    
## Fermeture de toutes les figures
plt.close('all')