
# coding: utf-8

# In[1]:

# Conditional Histograms for Boreal area of Variables conditioned to VCFF State
get_ipython().magic(u'matplotlib inline')
import matplotlib
#matplotlib.use('macosx')  # TkAgg Force mpl to use Tk backend
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import seaborn as sns
import datetime
import netCDF4
import sys
import math 
import matplotlib.patches as mpatches
# needs gridspec
from matplotlib import gridspec
from matplotlib.colors import ListedColormap
import palettable 
#sys.path.append("/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/") # Edited .zshrc profile with PYTHONPATH 
from pandas import DataFrame 
from collections import Counter
import pandas as pd



cmapbare = sns.light_palette("#332288", n_colors=5, reverse=False, as_cmap=True)
cmapsava = sns.light_palette("#44AA99", n_colors=5, reverse=False, as_cmap=True)
cmapsavab = sns.light_palette("#88CCEE", n_colors=5, reverse=False, as_cmap=True)
cmapfore = sns.light_palette("#117733", n_colors=5, reverse=False, as_cmap=True)

cmapfore_r = sns.light_palette("#117733", n_colors=5, reverse=True, as_cmap=True)
cmapbare_r = sns.light_palette("#332288", n_colors=5, reverse=True, as_cmap=True)
cmapsava_r = sns.light_palette("#44AA99", n_colors=5, reverse=True, as_cmap=True)
cmapsavab_r = sns.light_palette("#88CCEE", n_colors=5, reverse=True, as_cmap=True)

colbare = '#332288'
colsavab = '#88CCEE'
colsava = '#44AA99'
colfore = '#117733'

alt_colbare = '#999933'
alt_colsavab = '#CC6677' 
alt_colsavaf = '#AA4499'
alt_colfore = '#882255'


savabmap = ListedColormap(palettable.colorbrewer.diverging.PuOr_5_r.mpl_colors)
baremap = ListedColormap(palettable.colorbrewer.diverging.PuOr_5.mpl_colors)


def setcolor(x, color): # for basemap plots-> color parallels and meridians
    for m in x:
         for t in x[m][1]:
             t.set_color(color)







columns=["VCFF", "MAR", "Mean_Spring_SM", "Mean_Tmin", "PZI", "FF", "GDD0", "PTD", "ST", "State","Elev", "Lon", "Lat" ]
sns.set_style(style="white")


def altst_plot(m1, lonlat = np.zeros(8), dotsize = 20, alpha=0.8, marker = 'h'):
    #bare
    m1_sub_lons, m1_sub_lats = m1(lonlat[0],lonlat[1]) # convert lon lat to basemap coordinates
    cbs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=colbare,                 alpha=0.8, s=dotsize, zorder=3, label = 'Grassland / Open Woodland')
#m1.contourf(m1_sub_lons,m1_sub_lats, c_bare_0_bs,cmap=baremap, tri = True, alpha=1, s=100, zorder=6)

#savab
    m1_sub_lons, m1_sub_lats = m1(lonlat[2],lonlat[3])
    csb2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=colsavab,                 alpha=1, s=dotsize, zorder=4, label = "Open Woodland / Grassland")
#m1.contourf(m1_sub_lons,m1_sub_lats, c_sava_0_bs,cmap=savabmap, tri = True, alpha=0.8, s=50, zorder=6)

#savaf
    m1_sub_lons, m1_sub_lats = m1(lonlat[4],lonlat[5])
    csf2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=colsava,                 alpha=0.8, s=dotsize, zorder=4, label = "Open Woodland / Forest")


#fore
    m1_sub_lons, m1_sub_lats = m1(lonlat[6],lonlat[7])
    cfs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=colfore,                 alpha=0.8, s=dotsize, zorder=4, label = "Forest / Open Woodland")
    
    return cbs2, csb2, csf2, cfs2

def st_plot(m1, lonlat = np.zeros(8), dotsize = 20, alpha=0.8, marker = 'h'):
    #bare
    m1_sub_lons, m1_sub_lats = m1(lonlat[0],lonlat[1]) # convert lon lat to basemap coordinates
    cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=alt_colbare,                 alpha=0.8, s=dotsize, zorder=3, label = 'Grassland / Open Woodland')
#m1.contourf(m1_sub_lons,m1_sub_lats, c_bare_0_bs,cmap=baremap, tri = True, alpha=1, s=100, zorder=6)

#savab
    m1_sub_lons, m1_sub_lats = m1(lonlat[2],lonlat[3])
    csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=alt_colsavab,                 alpha=1, s=dotsize, zorder=4, label = "Open Woodland / Grassland")
#m1.contourf(m1_sub_lons,m1_sub_lats, c_sava_0_bs,cmap=savabmap, tri = True, alpha=0.8, s=50, zorder=6)

#savaf
    m1_sub_lons, m1_sub_lats = m1(lonlat[4],lonlat[5])
    csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=alt_colsavaf,                 alpha=0.8, s=dotsize, zorder=4, label = "Open Woodland / Forest")


#fore
    m1_sub_lons, m1_sub_lats = m1(lonlat[6],lonlat[7])
    cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker=marker, c=alt_colfore,                 alpha=0.8, s=dotsize, zorder=4, label = "Forest / Open Woodland")
    
    return cbs, csb, csf, cfs


na_e_dist_lonlat = np.load('na_e_dist_lolat.npy')
na_e_alt1_lonlat = np.load('na_e_alt1_lolat.npy')
na_e_alt2_lonlat = np.load('na_e_alt2_lolat.npy')
na_w_dist_lonlat = np.load('na_w_dist_lolat.npy')
na_w_alt1_lonlat = np.load('na_w_alt1_lolat.npy')
na_w_alt2_lonlat = np.load('na_w_alt2_lolat.npy')
ea_w_dist_lonlat = np.load('ea_w_dist_lolat.npy')
ea_w_alt1_lonlat = np.load('ea_w_alt1_lolat.npy')
ea_w_alt2_lonlat = np.load('ea_w_alt2_lolat.npy')
ea_e_dist_lonlat = np.load('ea_e_dist_lolat.npy')
ea_e_alt1_lonlat = np.load('ea_e_alt1_lolat.npy')
ea_e_alt2_lonlat = np.load('ea_e_alt2_lolat.npy')



# In[9]:

# possible states, but with different projection
fig = plt.figure(figsize=(18,10), facecolor='white')
ax = fig.add_subplot(111)
#ax.set_title("Possible Alternative Tree Cover States", fontsize=18)
m1 = Basemap(projection='npaeqd',boundinglat=45,lon_0=-90,resolution='l', round=True)

# map extra features # labels are at intersection the left, right, top or bottom of the plot
par = m1.drawparallels([45.,55.,65.,75.,90.],labels=[False,False,False,False], color='black', fontsize=24)
mer = m1.drawmeridians(np.arange(-180.,179.,45.),labels=[False,False,True,True], color='black', fontsize=24)

m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'
m1.drawcoastlines(linewidth=1.5)
m1.fillcontinents(color='white',lake_color='black')

# set the color of the parallel and meridian lines
setcolor(par,'black')
setcolor(mer,'black')



# fire disturbed
dotsize = 10
mar = 'h'
altst_plot(m1=m1, lonlat = na_e_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
altst_plot(m1=m1, lonlat = na_w_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
cbs2, csb2, csf2, cfs2 = altst_plot(m1=m1, lonlat = ea_e_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
altst_plot(m1=m1, lonlat = ea_w_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)

# nondist
st_plot(m1=m1, lonlat = na_e_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = na_w_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_e_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_w_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)


st_plot(m1=m1, lonlat = na_e_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = na_w_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
cbs, csb, csf, cfs = st_plot(m1=m1, lonlat = ea_e_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_w_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)


legend1 = plt.legend((cbs, csb, csf, cfs), ('Treeless / Open Woodland', 'Open Woodland / Treeless', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1.05, 0.3), markerscale = 5, ncol = 1, fontsize = 24, fancybox = True, frameon=True,                      title = "Fire Disturbed Tree Cover States")

legend2 = plt.legend((cbs2, csb2, csf2, cfs2), ('Treeless / Open Woodland', 'Open Woodland / Treeless', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1.05, 0.7), markerscale = 5, ncol = 1, fontsize = 24, fancybox = True, frameon=True,                      title = "Equivalent Tree Cover States")
legend2.get_title().set_fontsize('24') #legend 'Title' fontsize
legend1.get_title().set_fontsize('24') #legend 'Title' fontsize

plt.gca().add_artist(legend1)
plt.subplots_adjust(left=0.001)
#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Meeting_Victor/Alternative_States/New_Classification_with_SoilTexture_and_Disturbances/PanBoreal_Possible_Alternative_10.pdf', dpi=600, transparent = True)


plt.show()


# In[41]:

fig = plt.figure(figsize=(12,9), facecolor='white')

ax = fig.add_subplot(111)
#ax.set_title("Same Variables Tree Cover States over Eastern North America", fontsize=18)
#m = Basemap(projection='ortho',lon_0=-125,lat_0=60,resolution='l') #West
#m = Basemap(projection='ortho',lon_0=95.,lat_0=60,resolution='l') #East
#width = m.urcrnrx - m.llcrnrx
#height = m.urcrnry - m.llcrnry
# let's zoom in the plot
#coef = 0.35
#width = width*coef
#height = height*coef
#m1 = Basemap(projection='ortho',lon_0=75.,lat_0=60.,resolution='l',
 #           llcrnrx=-0.5*width,llcrnry=-0.5*height,urcrnrx=0.5*width,urcrnry=0.5*height)# same as previous but with new extension
# map extra features# labels are at intersection the left, right, top or bottom of the plot

#m1= Basemap(projection='cyl',llcrnrlat=45,urcrnrlat=80,\
#            llcrnrlon=70,urcrnrlon=135,resolution='l')


m1 = Basemap(width=6000000,height=4000000,
            resolution='l',projection='laea',\
            lat_ts=50,lat_0=65,lon_0=90.)

par = m1.drawparallels([45.,55.,65.,75.,90.],labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,179.,15.),labels=[False,False,False,True], color='black', fontsize=14)

m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'02c5e9
m1.drawcoastlines(linewidth=1.5)
m1.fillcontinents(color='white',lake_color='black')

# set the color of the parallel and meridian lines
setcolor(par,'black')
setcolor(mer,'black')



# fire disturbed
dotsize = 25
mar = 'h'
altst_plot(m1=m1, lonlat = na_e_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
altst_plot(m1=m1, lonlat = na_w_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
cbs2, csb2, csf2, cfs2 = altst_plot(m1=m1, lonlat = ea_e_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
altst_plot(m1=m1, lonlat = ea_w_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)

# nondist
st_plot(m1=m1, lonlat = na_e_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = na_w_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_e_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_w_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)


st_plot(m1=m1, lonlat = na_e_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = na_w_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
cbs, csb, csf, cfs = st_plot(m1=m1, lonlat = ea_e_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_w_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)










#legend1 = plt.legend((cbs, csb, csf, cfs), ('Treeless / Open Woodland', 'Open Woodland / Treeless', 'Open Woodland / Forest', 'Forest / Open Woodland') ,\
#                     loc='center left', bbox_to_anchor=(0.05, 0.3), markerscale = 4, ncol = 1, fontsize = 20, fancybox = True, frameon=True, \
#                     title = "Fire Disturbed Tree Cover States")

#legend2 = plt.legend((cbs2, csb2, csf2, cfs2), ('Treeless / Open Woodland', 'Open Woodland / Treeless', 'Open Woodland / Forest', 'Forest / Open Woodland') ,\
#                     loc='center left', bbox_to_anchor=(0.05, 0.7), markerscale = 4, ncol = 1, fontsize = 20, fancybox = True, frameon=True, \
#                     title = "Equivalent Tree Cover States")
#legend2.get_title().set_fontsize('20') #legend 'Title' fontsize
#legend1.get_title().set_fontsize('20') #legend 'Title' fontsize

#plt.gca().add_artist(legend1)
#plt.subplots_adjust(left=0.001)
#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/My_Papers/Copernicus_LaTeX_Package/pic/Eurasia_Possible_Alternative.pdf', dpi=600, transparent = False)


plt.show()


# In[76]:

fig = plt.figure(figsize=(16,9), facecolor='white')

ax = fig.add_subplot(111)


#ax.set_title("Same Variables Tree Cover States over Eastern North America", fontsize=18)

# map extra features# labels are at intersection the left, right, top or bottom of the plot

m1 = Basemap(width=6000000,height=4000000,
            resolution='l',projection='laea',\
            lat_ts=50,lat_0=65,lon_0=-95.)

par = m1.drawparallels([45.,55.,65.,75.,90.],labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-130.,179.,20.),labels=[False,False,False,True], color='black', fontsize=14)

m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'02c5e9
m1.drawcoastlines(linewidth=1.5)
m1.fillcontinents(color='white',lake_color='black')

# set the color of the parallel and meridian lines
setcolor(par,'black')
setcolor(mer,'black')



# fire disturbed
dotsize = 25
mar = 'h'
altst_plot(m1=m1, lonlat = na_e_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
altst_plot(m1=m1, lonlat = na_w_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
cbs2, csb2, csf2, cfs2 = altst_plot(m1=m1, lonlat = ea_e_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
altst_plot(m1=m1, lonlat = ea_w_dist_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)

# nondist
st_plot(m1=m1, lonlat = na_e_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = na_w_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_e_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_w_alt1_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)


st_plot(m1=m1, lonlat = na_e_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = na_w_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
cbs, csb, csf, cfs = st_plot(m1=m1, lonlat = ea_e_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)
st_plot(m1=m1, lonlat = ea_w_alt2_lonlat, dotsize = dotsize, alpha=0.8, marker = mar)


fig = plt.figure(figsize=(16,9), facecolor='white')

ax = fig.add_subplot(111)


legend1 = plt.legend((cbs, csb, csf, cfs), ('Treeless / Open Woodland', 'Open Woodland / Treeless', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='upper center', bbox_to_anchor=(0.5, 0.10), markerscale = 2, ncol = 4, fontsize = 12, fancybox = True, frameon=False,                      title = "Fire Disturbed Tree Cover States", borderaxespad=0.1)

legend2 = plt.legend((cbs2, csb2, csf2, cfs2), ('Treeless / Open Woodland', 'Open Woodland / Treeless', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='upper center', bbox_to_anchor=(0.5, 0.2), markerscale = 2, ncol = 4, fontsize = 12, fancybox = True, frameon=False,                      title = "Equivalent Tree Cover States")
legend2.get_title().set_fontsize('14') #legend 'Title' fontsize
legend1.get_title().set_fontsize('14') #legend 'Title' fontsize

plt.gca().add_artist(legend1)
#plt.gca().add_artist(legend2)

#plt.subplots_adjust(left=0.001)
#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/My_Papers/Copernicus_LaTeX_Package/pic/Legend_Possible_Alternative.pdf', dpi=600, transparent = False)


plt.show()


# In[ ]:



