
# 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 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)




# function to create the delimiters for the classification given the arrays of bare, savanna, forest
# it takes the second and second to last percentile and then divide interval in between them
# function to create the delimiters for the classification given the arrays of bare, savanna, forest
# it takes the second and second to last percentile and then divide interval in between them
def delimiters_generator(bare, savanna, forest, indexes = [1,2,3,4,5], intervals = 6):
    percentiles = [] # list that will contain all the 10th and 90th percentiles related to the variables from indexes 
    for i in indexes:
        percentiles_dummy = []
        percentiles_dummy.append(np.percentile(forest[:,i],90))
        percentiles_dummy.append(np.percentile(forest[:,i],10))
        percentiles_dummy.append(np.percentile(savanna[:,i],90))
        percentiles_dummy.append(np.percentile(savanna[:,i],10))
        percentiles_dummy.append(np.percentile(bare[:,i],90))
        percentiles_dummy.append(np.percentile(bare[:,i],10))
        percentiles_dummy = sorted(percentiles_dummy, key=float) # sort the list in ascending order
        percentiles.append(percentiles_dummy)
    delimiters = [np.linspace(percentiles[i][1], percentiles[i][4],intervals-1) for i in [0,1,2]]
    for j in [3,4]:
        first = np.nonzero(percentiles[j])[0][0]
        second = np.nonzero(percentiles[j])[0][-2]
        if (j==3): delimiters.append(np.array([percentiles[j][first]/2, percentiles[j][first], percentiles[j][second]]))
        else: delimiters.append(np.array([0.,percentiles[j][first]/2, percentiles[j][first], percentiles[j][second]]))

    return delimiters

# the classification algorithm, looks at the delimiters and creates bins for the input variables
def ea_east_classifier(variables_array, delimiters):
    classified = np.zeros(shape = variables_array[:,0:3].shape)

    
    for i in range(classified.shape[0]):
        classified[i,1]=variables_array[i,-2] # lon
        classified[i,2]=variables_array[i,-1] # lat
        # using mean minimum temperature
        if (variables_array[i,1]<=delimiters[0][0]): 
            classified[i,0]=1
        elif (variables_array[i,1]<=delimiters[0][1]):
            classified[i,0]=2
        elif (variables_array[i,1]<=delimiters[0][2]):
            classified[i,0]=3
        elif (variables_array[i,1]<=delimiters[0][3]):
            classified[i,0]=4
        elif (variables_array[i,1]<=delimiters[0][4]):
            classified[i,0]=5
        elif (variables_array[i,1]>delimiters[0][4]):
            classified[i,0]=6
        
        # using mean spring soil moisture
        if (variables_array[i,2]<=delimiters[1][0]):
            classified[i,0]+=0
        elif (variables_array[i,2]<=delimiters[1][1]):
            classified[i,0]+=6
        elif (variables_array[i,2]<=delimiters[1][2]):
            classified[i,0]+=12
        elif (variables_array[i,2]<=delimiters[1][3]):
            classified[i,0]+=18
        elif (variables_array[i,2]<=delimiters[1][4]):
            classified[i,0]+=24
        elif (variables_array[i,2]>delimiters[1][4]):
            classified[i,0]+=30
        
        # using mean annual rainfall
        if (variables_array[i,3]<=delimiters[2][0]):
            classified[i,0]+=0
        elif (variables_array[i,3]<=delimiters[2][1]):
            classified[i,0]+=36
        elif (variables_array[i,3]<=delimiters[2][2]):
            classified[i,0]+=72
        elif (variables_array[i,3]<=delimiters[2][3]):
            classified[i,0]+=108
        elif (variables_array[i,3]<=delimiters[2][4]):
            classified[i,0]+=144
        elif (variables_array[i,3]>delimiters[2][4]):
            classified[i,0]+=180

           
        # using permafrost zonation index
        if (variables_array[i,4]<=delimiters[3][0]):
            classified[i,0]+=0
        elif (variables_array[i,4]<=delimiters[3][1]):
            classified[i,0]+=216
        elif (variables_array[i,4]<=delimiters[3][2]):
            classified[i,0]+=432
        elif (variables_array[i,4]>delimiters[3][2]):
            classified[i,0]+=648

        
        # using soil texture
        # better put soil texture before fire, otherwise it will be impossible to check the fire as disturing agent
        if (np.any(variables_array[i,8]==[1,2,12])):
            classified[i,0]+=0. 
        elif (np.any(variables_array[i,8]==[3,4,5,6,7,8])):
            classified[i,0]+=1000. 
        elif (np.any(variables_array[i,8]==[9,10,11,13,14])):
            classified[i,0]+=2000.
        else:
            classified[i,0]+=3000. 

        
        # using fire frequency
        if (variables_array[i,5]<=delimiters[4][0]):
            classified[i,0]+=0
        elif (variables_array[i,5]<=delimiters[4][1]):
            classified[i,0]+=3864
        elif (variables_array[i,5]<=delimiters[4][2]):
            classified[i,0]+=3864*2
        elif (variables_array[i,5]<=delimiters[4][3]):
            classified[i,0]+=3864*3
        elif (variables_array[i,5]>delimiters[4][3]):
            classified[i,0]+=3864*4
 
            
            
    return classified




def generate_legend(fig,c_bare_2_bs, c_sava_2_bs, c_fore_2_fs, c_sava_2_fs, save_2_bs, save_2_sb, save_2_fs, save_2_sf, cmaps = [ListedColormap(palettable.colorbrewer.diverging.PuOr_5.mpl_colors), ListedColormap(palettable.colorbrewer.diverging.PuOr_5_r.mpl_colors), 'summer', 'winter'] ):
    un = [np.unique(c_bare_2_bs), np.unique(c_sava_2_bs), np.unique(c_fore_2_fs), np.unique(c_sava_2_fs)]

    #cmaps = [baremap, savabmap, 'summer', 'winter'] #['Reds_r', 'spring', 'summer', 'winter']
    rep = [[np.random.rand(1)[0] for elem in un[i]] for i in range(4)]
    rop = [[np.random.rand(1)[0] for elem in un[i]] for i in range(4)]
    ck = [plt.scatter(rep[i], rop[i], c=un[i], cmap=cmaps[i], s=500, label=rep[i]) for i in range(4)]#'terrain_r'

#patches = [mpatches.Patch(facecolor=plt.get(ck[2], 'facecolors')[i], label=labels[i]) for i in range(len(plt.get(ck[2], 'facecolors')))]
#patch = [mpatches.Patch(facecolor=plt.get(ck[2], 'facecolors')[0], label = 'test')]
#plt.legend(handles=patch, loc='lower left')
    plt.draw()
    ax1 = fig.add_subplot(111)

    col = []
    for j in range(4):
        col.append([plt.get(ck[j], 'facecolors')[i] for i in range(len(plt.get(ck[j], 'facecolors')))])
    col
    len_tot = 0
    for elem in col: 
        len_tot += len(elem)
    len_tot
    x = np.linspace(0, 1, 10)
    patches = []
    ind = 0
#temporary labels, fix them later
    s_1a = [np.unique(save_2_bs)[i] for i in range(len(np.unique(save_2_bs))) if (np.in1d(np.unique(save_2_bs),np.unique(c_bare_2_bs))[i])]
    s_2a = [np.unique(save_2_fs)[i] for i in range(len(np.unique(save_2_fs))) if (np.in1d(np.unique(save_2_fs),np.unique(c_fore_2_fs))[i])]
    
    s_1 = [np.unique(save_2_sb)[i] for i in range(len(np.unique(save_2_sb))) if (np.in1d(np.unique(save_2_sb),np.unique(c_sava_2_bs))[i])]
    s_2 = [np.unique(save_2_sf)[i] for i in range(len(np.unique(save_2_sf))) if (np.in1d(np.unique(save_2_sf),np.unique(c_sava_2_fs))[i])]
    lab = [str(i) for i in np.concatenate((s_1a, s_1, s_2a, s_2), axis=0)]
    labels = []
    for i in range(len(s_1a)):
        labels.append('Grassland/Open Woodland '+lab[i])
    for i in range(len(s_1)):
        labels.append('Open Woodland/Grassland '+lab[i+len(s_1a)])
    for i in range(len(s_2a)):
        labels.append('Forest/Open Woodland '+lab[i+len(s_1a)+len(s_1)])
    for i in range(len(s_2)):
        labels.append('Open Woodland/Forest '+lab[i+len(s_1a)+len(s_1)+len(s_2a)])
    for j in range(len(col)):
        for i in range(len(col[j])):
            plt.plot(x, np.random.rand(1) * x + np.random.rand(1), color = col[j][i], label='$y = {i}x + {i}$'.format(i=i))
            patches.append(mpatches.Patch(facecolor=col[j][i], label=labels[ind]))
            ind +=1
    plt.legend(handles = patches, loc='best')
    plt.draw()
    return patches

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)


from Boreal_Variables import bare_ea_east as bare
from Boreal_Variables import savanna_ea_east as savanna
from Boreal_Variables import forest_ea_east as forest 

from Boreal_Variables import bare_ea_east as bare_tea_e
from Boreal_Variables import savanna_ea_east as savanna_tea_e
from Boreal_Variables import forest_ea_east as forest_tea_e

from Boreal_Variables import bare_ea_west_russia as bare_tea_wr
from Boreal_Variables import savanna_ea_west_russia as savanna_tea_wr
from Boreal_Variables import forest_ea_west_russia as forest_tea_wr

from Boreal_Variables import bare_ea_europe as bare_tea_eu
from Boreal_Variables import savanna_ea_europe as savanna_tea_eu
from Boreal_Variables import forest_ea_europe as forest_tea_eu

from Boreal_Variables import bare_na_east_update as bare_tna_e
from Boreal_Variables import savanna_na_east_update as savanna_tna_e
from Boreal_Variables import forest_na_east_update as forest_tna_e

from Boreal_Variables import bare_na_west_update as bare_tna_w
from Boreal_Variables import savanna_na_west_update as savanna_tna_w
from Boreal_Variables import forest_na_west_update as forest_tna_w






dataforr_ea_east = np.loadtxt("/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/boreal_ea_east_database.txt") # columns=["VCFF", "MAR", "MASM", "Tmin", "PZI", "FF", "GDD0", "State", "Elev", "Lon", "Lat"]
columns=["VCFF", "MAR", "Mean_Spring_SM", "Mean_Tmin", "PZI", "FF", "GDD0", "PTD", "ST", "State","Elev", "Lon", "Lat" ]
dataframe_ea_east = DataFrame(dataforr_ea_east[:,0:13], columns=columns[0:13])


sns.set_style(style="white")

beat = len(bare_tea_e)+len(bare_tea_wr)+len(bare_tea_eu)
seat = len(savanna_tea_e)+len(savanna_tea_eu)+len(savanna_tea_wr)
feat = len(forest_tea_e)+len(forest_tea_eu)+len(forest_tea_wr)


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='h', c=alt_colbare,                 alpha=0.8, s=20, 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='h', c=alt_colsavab,                 alpha=1, s=20, 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='h', c=alt_colsavaf,                 alpha=0.8, s=20, 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='h', c=alt_colfore,                 alpha=0.8, s=20, zorder=4, label = "Forest / Open Woodland")
    
    return
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
    cbs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=20, 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='h', c=colsavab,                 alpha=1, s=20, 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='h', c=colsava,                 alpha=0.8, s=20, 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='h', c=colfore,                 alpha=0.8, s=20, zorder=4, label = "Forest / Open Woodland")
    
    return

#lonlat = [lon_distubare2,lat_distubare2,lon_distusavab2,lat_distusavab2, lon_distusavaf2,lat_distusavaf2, lon_distufore2,lat_distufore2]
#lonlat2 = [lon_bare_1_sb,lat_bare_1_sb, lon_sava_1_sb,lat_sava_1_sb, lon_sava_1_sf,lat_sava_1_sf, lon_fore_1_sf,lat_fore_1_sf]
#lonlat3 = [lon_bare_2_bs,lat_bare_2_bs, lon_sava_2_bs,lat_sava_2_bs, lon_sava_2_fs,lat_sava_2_fs, lon_fore_2_fs,lat_fore_2_fs]
#np.save('ea_e_dist_lolat', lonlat)
#np.save('ea_e_alt1_lolat', lonlat2)
#np.save('ea_e_alt2_lolat', lonlat3)
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[15]:

# set delimiters and run classification # if you load the pickles, don't run this

# delimiters west_russia ea (18)
#delimiters = delimiters_generator(bare, savanna, forest)

#col_class = ["state", "lon", "lat"]
#bare_classified = DataFrame(ea_east_classifier(bare, delimiters), columns = col_class)
#savanna_classified = DataFrame(ea_east_classifier(savanna, delimiters), columns = col_class)
#forest_classified = DataFrame(ea_east_classifier(forest, delimiters), columns = col_class)

#bare_classified.to_pickle('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/bare_ea_east_classified_with_SoilT.pkl')
#savanna_classified.to_pickle('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/savanna_ea_east_classified_with_SoilT.pkl')
#forest_classified.to_pickle('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/forest_ea_east_classified_with_SoilT.pkl')

# delimiters


# In[2]:

delimiters = delimiters_generator(bare, savanna, forest)
for elem in delimiters: print elem


# In[20]:

# load dataframes of classified
bare_classified = pd.read_pickle('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/bare_ea_east_classified_with_SoilT.pkl')
savanna_classified = pd.read_pickle('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/savanna_ea_east_classified_with_SoilT.pkl')
forest_classified = pd.read_pickle('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/forest_ea_east_classified_with_SoilT.pkl')

# count the intersecting values with at least a minimum amount of repetitions
sb0count, sbcount = [], [] # counters for savanna-bare
sf0count, sfcount = [], [] # counters for savanna-forest
b0count, bcount = [], [] # counters for bare
f0count, fcount = [], [] # counters for forest

index_0_bs, save_0_bs = [], [] # index to eliminate or keep of elements which are in the intersections and appear a minimum amount of times
index_1_sb, save_1_sb = [], [] # for all variables are the same (0), all variables are the same except fire for forest and bare (2)
index_2_bs, save_2_bs = [], [] # all the variables are the same except fire for savanna (1)   
index_0_fs, save_0_fs = [], [] # [sf means savanna fixed, forest disturbed and vice versa, except in the 0 case - and the intersection is between forest and savanna]
index_1_sf, save_1_sf = [], []
index_2_fs, save_2_fs = [], []

minvalue_b, minvalue_s, minvalue_f = int(0.005*len(bare_tea_e)), int(0.005*len(savanna_tea_e)), int(0.01*len(forest_tea_e)) # min values of repetitions for states
print minvalue_b, minvalue_s, minvalue_f
#minvalue_b, minvalue_s, minvalue_f = 3,3,3
fire_subtracter = 3864 # the value for fire in the classification 
template2 = "{0:18}|{1:7}|{2:12}|{3:5}|{4:12}|{5:13}" # template for printing function 
for st_number in np.unique(bare_classified["state"]): # loop through unique elements of bare_classified
    b0count.append(Counter(bare_classified["state"])[st_number]) # count repetitions of each element in bare
    sb0count.append(Counter(savanna_classified["state"])[st_number]) # and in savanna
    bc, sb = [], []
    for subtract_mult_bare in range(1,4): # loop in the possible bins for fire
        bc.append(Counter(bare_classified["state"]-fire_subtracter*subtract_mult_bare)[st_number]) # count repetitions
        sb.append(Counter(savanna_classified["state"]-fire_subtracter*subtract_mult_bare)[st_number])
    bcount.append(bc)
    sbcount.append(sb)
        #print bcount #l'intersezione deve anche essere lista

        
        
for i in range(len(b0count)): # loop number of through counted elements
    if (b0count[i]>=10*sb0count[i] or sb0count[i]>=10*b0count[i]): index_0_bs.append(i) # if one counter is at least 10folds the other, discard
    elif (b0count[i]>minvalue_b and sb0count[i]>minvalue_s): # if both counters are above thresholds, save the corresponding classified state 
        print template2.format("bare/op wood state", np.unique(bare_classified["state"])[i], "b count", b0count[i], "op wo count", sb0count[i])
        save_0_bs.append(np.unique(bare_classified["state"])[i]) # the saved classified states in the case where all variables are the same
    else: index_0_bs.append(i)

for i in range(len(b0count)): # same as previous loop
    if (sb0count[i]>minvalue_s and (np.any(np.array(bcount[i])>minvalue_b))): # counts if the standard savanna counter and any one of the disturbed ones are > thresholds
        print template2.format("op wood state", np.unique(bare_classified["state"])[i], "op wo count", sb0count[i], "b count", bcount[i])
        save_1_sb.append(np.unique(bare_classified["state"])[i]) # save the savab fixed and bare disturbed states 
    else: index_1_sb.append(i)
        
for i in range(len(b0count)): # same as previous, but bare is fixed and savanna is fire disturbed
    if (np.all(b0count[i]>=10*np.array(sbcount[i])) or np.any(np.array(sbcount[i])>=10*(b0count[i]+1))): index_2_bs.append(i)
    elif (b0count[i]>minvalue_b and (np.any(np.array(sbcount[i])>minvalue_s))):
        print template2.format("bare state", np.unique(bare_classified["state"])[i], "b count", b0count[i], "op wo count", sbcount[i])
        save_2_bs.append(np.unique(bare_classified["state"])[i]) # bare fixed and savab disturbed
    else: index_2_bs.append(i)
        
print " "          
        
# same as the previous loops, but this time between savanna and forest states     
for st_number in np.unique(forest_classified["state"]):
    f0count.append(Counter(forest_classified["state"])[st_number])
    sf0count.append(Counter(savanna_classified["state"])[st_number])
    fc, sf = [], []
    for subtract_mult_fore in range(1,4):
        fc.append(Counter(forest_classified["state"]-fire_subtracter*subtract_mult_fore)[st_number])
        sf.append(Counter(savanna_classified["state"]-fire_subtracter*subtract_mult_fore)[st_number])
    fcount.append(fc)
    sfcount.append(sf)
    
f0total = 0     
sf0total = 0
for i in range(len(f0count)):
    if (f0count[i]>=10*sf0count[i] or sf0count[i]>=10*f0count[i]): index_0_fs.append(i)
    elif (f0count[i]>minvalue_f and sf0count[i]>minvalue_s):
        print template2.format("fore/op wood state", np.unique(forest_classified["state"])[i], "f count", f0count[i], "op wo count", sf0count[i])
        save_0_fs.append(np.unique(forest_classified["state"])[i])
        f0total += f0count[i]
        sf0total += sf0count[i]
    else: index_0_fs.append(i)

        
ftotal = 0
sftotal = 0
for i in range(len(f0count)):
      
    #if (np.any(np.array(bcount[i])>=10*sb0count[i]) or np.any(sb0count[i]>=10*np.array(b0count[i]))): index_1_sf.append(i)
    if (sf0count[i]>minvalue_s and (np.any(np.array(fcount[i])>minvalue_f))):
        print template2.format("op wood state", np.unique(forest_classified["state"])[i], "op wo count", sf0count[i], "f count", fcount[i])
        save_1_sf.append(np.unique(forest_classified["state"])[i])
        ftotal += np.sum(fcount[i])
        sftotal += sf0count[i]
    else: index_1_sf.append(i)

f2total = 0 
sf2total = 0
for i in range(len(f0count)):
      
    #if (bcount[i]>=10*sb0count[i] or sb0count[i]>=10*b0count[i]): index_0.append(i)
    if (f0count[i]>minvalue_f and (np.any(np.array(sfcount[i])>minvalue_s))):
        print template2.format("forest state", np.unique(forest_classified["state"])[i], "f count", f0count[i], "op wo count", sfcount[i])
        save_2_fs.append(np.unique(forest_classified["state"])[i])
        sf2total += np.sum(sfcount[i])
        f2total += f0count[i]
    else: index_2_fs.append(i)
    

    
## here the indexes contain the positions of elements in the unique vectors of the classified states  that should be elimited in the mask
# to get the actual element you should check if the element corresponds to a certain state which corresponds to a certain ordered element from the
# indexes 
# so to do so you need to save the state that you don't want to erase (because it's faster this way) and then
# create the masks by looking at which elements of the original _classified are also present in the filtered save_ 

# print results: case - state - type of counter - count - type of counter - count


# In[21]:

sums = f0total, sf0total, sftotal, sf2total, ftotal, f2total
total = np.sum(sums)
print sums, total, len(forest), len(savanna)


# In[22]:

# masks for where all the datasets are used for classification and I cannot separate certains states
# let's create the masks : 4 level_0 masks (1 bare, 1 savab, 1 savaf, 1 fore) and corresponding lonlat masks

mask_bare_0_bs = np.in1d(bare_classified['state'],save_0_bs) # create mask of intersection
mask_sava_0_bs = np.in1d(savanna_classified['state'],save_0_bs) # keeps the element of classified that are in the save vector
mask_sava_0_fs = np.in1d(savanna_classified['state'],save_0_fs)
mask_fore_0_fs = np.in1d(forest_classified['state'],save_0_fs)


c_bare_0_bs = np.ma.masked_array(bare_classified['state'], ~mask_bare_0_bs) #apply inverse mask
lon_bare_0_bs = np.ma.masked_array(bare_classified['lon'], ~mask_bare_0_bs) # creates a masked array same as the classified but according to the mask
lat_bare_0_bs = np.ma.masked_array(bare_classified['lat'], ~mask_bare_0_bs)

c_sava_0_bs = np.ma.masked_array(savanna_classified['state'], ~mask_sava_0_bs) #apply inverse mask
lon_sava_0_bs = np.ma.masked_array(savanna_classified['lon'], ~mask_sava_0_bs)
lat_sava_0_bs = np.ma.masked_array(savanna_classified['lat'], ~mask_sava_0_bs)

c_fore_0_fs = np.ma.masked_array(forest_classified['state'], ~mask_fore_0_fs) #apply inverse mask
lon_fore_0_fs = np.ma.masked_array(forest_classified['lon'], ~mask_fore_0_fs)
lat_fore_0_fs = np.ma.masked_array(forest_classified['lat'], ~mask_fore_0_fs)

c_sava_0_fs = np.ma.masked_array(savanna_classified['state'], ~mask_sava_0_fs) #apply inverse mask
lon_sava_0_fs = np.ma.masked_array(savanna_classified['lon'], ~mask_sava_0_fs)
lat_sava_0_fs = np.ma.masked_array(savanna_classified['lat'], ~mask_sava_0_fs)


# masks for fire as disturbance
# let's create the 8 masks : 
# 1 for bare fixed and 1 bare fire disturbed
# 1 for savab fixed and 1 savab disturbed
# 1 for savaf fixed and 1 savaf disturbed
# 1 for fore fixed and 1 fore disturbed 

# create mask of intersection - bare fixed
mask_bare_2_bs = np.in1d(bare_classified['state'],save_2_bs)  # 2 = bare fixed for bare
save_2_sb = np.ravel([elem+fire_subtracter*np.array((1,2,3)) for elem in save_2_bs]) # consider all elements modulo firesubtracter
mask_sava_2_bs = np.in1d(savanna_classified['state'],save_2_sb)  # 2 = bare fixed for sava
# create the inverse mask and the mask for lon lat from the one just created - bare
c_bare_2_bs = np.ma.masked_array(bare_classified['state'], ~mask_bare_2_bs) #apply inverse mask
lon_bare_2_bs = np.ma.masked_array(bare_classified['lon'], ~mask_bare_2_bs)
lat_bare_2_bs = np.ma.masked_array(bare_classified['lat'], ~mask_bare_2_bs)
# create the inverse mask and the mask for lon lat from the one just created - savanna
c_sava_2_bs = np.ma.masked_array(savanna_classified['state'], ~mask_sava_2_bs) #apply inverse mask
lon_sava_2_bs = np.ma.masked_array(savanna_classified['lon'], ~mask_sava_2_bs)
lat_sava_2_bs = np.ma.masked_array(savanna_classified['lat'], ~mask_sava_2_bs)

# create mask of intersection - fore fixed 
mask_fore_2_fs = np.in1d(forest_classified['state'],save_2_fs)  # 2 = fore fixed
save_2_sf = np.ravel([elem+fire_subtracter*np.array((1,2,3)) for elem in save_2_fs])
mask_sava_2_fs = np.in1d(savanna_classified['state'],save_2_sf)  # 2 = fore fixed
# create the inverse mask and the mask for lon lat from the one just created - forest
c_fore_2_fs = np.ma.masked_array(forest_classified['state'], ~mask_fore_2_fs) #apply inverse mask
lon_fore_2_fs = np.ma.masked_array(forest_classified['lon'], ~mask_fore_2_fs)
lat_fore_2_fs = np.ma.masked_array(forest_classified['lat'], ~mask_fore_2_fs)
# create the inverse mask and the mask for lon lat from the one just created - savanna
c_sava_2_fs = np.ma.masked_array(savanna_classified['state'], ~mask_sava_2_fs) #apply inverse mask
lon_sava_2_fs = np.ma.masked_array(savanna_classified['lon'], ~mask_sava_2_fs)
lat_sava_2_fs = np.ma.masked_array(savanna_classified['lat'], ~mask_sava_2_fs)

# create mask of intersection - savab fixed
mask_sava_1_sb = np.in1d(savanna_classified['state'],save_1_sb)  # 1 = sava fixed
save_1_bs = np.ravel([elem+fire_subtracter*np.array((1,2,3)) for elem in save_1_sb])
mask_bare_1_sb = np.in1d(bare_classified['state'],save_1_bs)  # 1 = sava fixed
# create the inverse mask and the mask for lon lat from the one just created - bare
c_bare_1_sb = np.ma.masked_array(bare_classified['state'], ~mask_bare_1_sb) #apply inverse mask
lon_bare_1_sb = np.ma.masked_array(bare_classified['lon'], ~mask_bare_1_sb)
lat_bare_1_sb = np.ma.masked_array(bare_classified['lat'], ~mask_bare_1_sb)
# create the inverse mask and the mask for lon lat from the one just created - savanna
c_sava_1_sb = np.ma.masked_array(savanna_classified['state'], ~mask_sava_1_sb) #apply inverse mask
lon_sava_1_sb = np.ma.masked_array(savanna_classified['lon'], ~mask_sava_1_sb)
lat_sava_1_sb = np.ma.masked_array(savanna_classified['lat'], ~mask_sava_1_sb)

# create mask of intersection - savaf fixed
mask_sava_1_sf = np.in1d(savanna_classified['state'],save_1_sf)  # 1 = sava fixed
save_1_fs = np.ravel([elem+fire_subtracter*np.array((1,2,3)) for elem in save_1_sf])
mask_fore_1_sf = np.in1d(forest_classified['state'],save_1_fs)  # 1 = sava fixed
# create the inverse mask and the mask for lon lat from the one just created - forest
c_fore_1_sf = np.ma.masked_array(forest_classified['state'], ~mask_fore_1_sf) #apply inverse mask
lon_fore_1_sf = np.ma.masked_array(forest_classified['lon'], ~mask_fore_1_sf)
lat_fore_1_sf = np.ma.masked_array(forest_classified['lat'], ~mask_fore_1_sf)
# create the inverse mask and the mask for lon lat from the one just created - savanna
c_sava_1_sf = np.ma.masked_array(savanna_classified['state'], ~mask_sava_1_sf) #apply inverse mask
lon_sava_1_sf = np.ma.masked_array(savanna_classified['lon'], ~mask_sava_1_sf)
lat_sava_1_sf = np.ma.masked_array(savanna_classified['lat'], ~mask_sava_1_sf)

# create masks to disentangle the effects (remove duplicates)
mask_distufore = np.in1d(c_fore_0_fs,save_2_fs)
c_distufore = np.ma.masked_array(c_fore_0_fs, mask_distufore) #apply mask
lon_distufore = np.ma.masked_array(lon_fore_0_fs, mask_distufore)
lat_distufore = np.ma.masked_array(lat_fore_0_fs, mask_distufore)

mask_distusavaf = np.in1d(c_sava_0_fs,save_2_sf)
c_distusavaf = np.ma.masked_array(c_sava_0_fs, mask_distusavaf) #apply mask
lon_distusavaf = np.ma.masked_array(lon_sava_0_fs, mask_distusavaf)
lat_distusavaf = np.ma.masked_array(lat_sava_0_fs, mask_distusavaf)

mask_distufore2 = np.in1d(c_distufore,save_1_fs)
c_distufore2 = np.ma.masked_array(c_distufore, mask_distufore2) #apply mask
lon_distufore2 = np.ma.masked_array(lon_distufore, mask_distufore2)
lat_distufore2 = np.ma.masked_array(lat_distufore, mask_distufore2)

mask_distusavaf2 = np.in1d(c_distusavaf,save_1_sf)
c_distusavaf2 = np.ma.masked_array(c_distusavaf, mask_distusavaf2) #apply mask
lon_distusavaf2 = np.ma.masked_array(lon_distusavaf, mask_distusavaf2)
lat_distusavaf2 = np.ma.masked_array(lat_distusavaf, mask_distusavaf2)


mask_distubare = np.in1d(c_bare_0_bs,save_1_bs)
c_distubare = np.ma.masked_array(c_bare_0_bs, mask_distubare) #apply mask
lon_distubare = np.ma.masked_array(lon_bare_0_bs, mask_distubare)
lat_distubare = np.ma.masked_array(lat_bare_0_bs, mask_distubare)

mask_distusavab = np.in1d(c_sava_0_bs,save_1_sb)
c_distusavab = np.ma.masked_array(c_sava_0_bs, mask_distusavab) #apply mask
lon_distusavab = np.ma.masked_array(lon_sava_0_bs, mask_distusavab)
lat_distusavab = np.ma.masked_array(lat_sava_0_bs, mask_distusavab)

mask_distubare2 = np.in1d(c_distubare,save_2_bs)
c_distubare2 = np.ma.masked_array(c_distubare, mask_distubare2) #apply mask
lon_distubare2 = np.ma.masked_array(lon_distubare, mask_distubare2)
lat_distubare2 = np.ma.masked_array(lat_distubare, mask_distubare2)

mask_distusavab2 = np.in1d(c_distusavab,save_2_sb)
c_distusavab2 = np.ma.masked_array(c_distusavab, mask_distusavab2) #apply mask
lon_distusavab2 = np.ma.masked_array(lon_distusavab, mask_distusavab2)
lat_distusavab2 = np.ma.masked_array(lat_distusavab, mask_distusavab2)
# need to mask also the lon lat arrays



# In[23]:

fig = plt.figure(figsize=(12,12), facecolor='white')
#patches = generate_legend(fig,c_bare_1_sb, c_sava_1_sb, c_fore_2_fs, c_sava_2_fs, save_1_bs,save_1_sb,  save_2_fs, save_2_sf, cmaps=[cmapbare, cmapsava, cmapfore, 'winter'])

ax = fig.add_subplot(111)
ax.set_title("Fire Disturbed Tree Cover States over Northern East Eurasia", fontsize=18)
#m = Basemap(projection='ortho',lon_0=-125,lat_0=60,resolution='l') #West
m = Basemap(projection='ortho',lon_0=100.,lat_0=60,resolution='l') #East
width = m.urcrnrx - m.llcrnrx
height = m.urcrnry - m.llcrnry
# let's zoom in the plot, cause the original basemap is too zoomed out
coef = 0.45 # zooming factor (percentage of original dimensions)
width = width*coef
height = height*coef
m1 = Basemap(projection='ortho',lon_0=100.,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
par = m1.drawparallels(np.arange(-80.,81.,20.),labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,181.,20.),labels=[False,False,False,True], color='black', fontsize=14)
m1.drawmapboundary(fill_color='black')
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')

#bare1
m1_sub_lons, m1_sub_lats = m1(lon_bare_1_sb,lat_bare_1_sb)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=50, zorder=3, label = 'Grassland / Open Woodland')
#bare2
m1_sub_lons, m1_sub_lats = m1(lon_bare_2_bs,lat_bare_2_bs)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=50, zorder=3)
#savab1
m1_sub_lons, m1_sub_lats = m1(lon_sava_1_sb,lat_sava_1_sb)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=0.8, s=50, zorder=4, label = 'Open Woodland / Grassland')
#savab2
m1_sub_lons, m1_sub_lats = m1(lon_sava_2_bs,lat_sava_2_bs)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=0.8, s=50, zorder=4)
#savaf1
m1_sub_lons, m1_sub_lats = m1(lon_sava_1_sf,lat_sava_1_sf)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=50, zorder=4, label = 'Open Woodland / Forest')
#savaf2
m1_sub_lons, m1_sub_lats = m1(lon_sava_2_fs,lat_sava_2_fs)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=50, zorder=4)
#fore1
m1_sub_lons, m1_sub_lats = m1(lon_fore_1_sf,lat_fore_1_sf)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=50, zorder=4, label = 'Forest / Open Woodland')
#fore2
m1_sub_lons, m1_sub_lats = m1(lon_fore_2_fs,lat_fore_2_fs)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=50, zorder=4)

plt.legend(bbox_to_anchor=(0.5, -0.03), ncol = 2, loc='upper center', markerscale = 2, fontsize = 14, fancybox = True, frameon=True)
#legend = plt.legend(handles=patches, loc='upper center', bbox_to_anchor=(0.5, -0.03), ncol = 3, fontsize = 14, fancybox = True, frameon=True)


#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Meeting_Victor/Alternative_States/New_Classification_with_SoilTexture_and_Disturbances/Nea_East_Fire_Disturbed_States_002.pdf', dpi=600)


plt.show()


# In[24]:

fig = plt.figure(figsize=(15,10), facecolor='white')
#patches = generate_legend(fig,c_bare_1_sb, c_sava_1_sb, c_fore_2_fs, c_sava_2_fs, save_1_bs,save_1_sb,  save_2_fs, save_2_sf, cmaps=[cmapbare, cmapsava, cmapfore, 'winter'])

ax = fig.add_subplot(111)
ax.set_title("Possible Alternative Tree Cover States over Northern East Eurasia", fontsize=18)
#m = Basemap(projection='ortho',lon_0=-125,lat_0=60,resolution='l') #West
m = Basemap(projection='ortho',lon_0=100.,lat_0=60,resolution='l') #East
width = m.urcrnrx - m.llcrnrx
height = m.urcrnry - m.llcrnry
# let's zoom in the plot, cause the original basemap is too zoomed out
coef = 0.45 # zooming factor (percentage of original dimensions)
width = width*coef
height = height*coef
m1 = Basemap(projection='ortho',lon_0=100.,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
par = m1.drawparallels(np.arange(-80.,81.,20.),labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,181.,20.),labels=[False,False,False,True], color='black', fontsize=14)
m1.drawmapboundary(fill_color='black')
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')

#bare1
m1_sub_lons, m1_sub_lats = m1(lon_bare_1_sb,lat_bare_1_sb)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=50, zorder=3, label = 'Grassland / Open Woodland')
#bare2
m1_sub_lons, m1_sub_lats = m1(lon_bare_2_bs,lat_bare_2_bs)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=50, zorder=3)
#savab1
m1_sub_lons, m1_sub_lats = m1(lon_sava_1_sb,lat_sava_1_sb)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=0.8, s=50, zorder=4, label = 'Open Woodland / Grassland')
#savab2
m1_sub_lons, m1_sub_lats = m1(lon_sava_2_bs,lat_sava_2_bs)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=0.8, s=50, zorder=4)
#savaf1
m1_sub_lons, m1_sub_lats = m1(lon_sava_1_sf,lat_sava_1_sf)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=50, zorder=4, label = 'Open Woodland / Forest')
#savaf2
m1_sub_lons, m1_sub_lats = m1(lon_sava_2_fs,lat_sava_2_fs)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=50, zorder=4)
#fore1
m1_sub_lons, m1_sub_lats = m1(lon_fore_1_sf,lat_fore_1_sf)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=50, zorder=4, label = 'Forest / Open Woodland')
#fore2
m1_sub_lons, m1_sub_lats = m1(lon_fore_2_fs,lat_fore_2_fs)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=50, zorder=4)


# not fire disturbed
#bare
m1_sub_lons, m1_sub_lats = m1(lon_distubare2,lat_distubare2) # convert lon lat to basemap coordinates
cbs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colbare,                 alpha=0.8, s=50, 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(lon_distusavab2,lat_distusavab2)
csb2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colsavab,                 alpha=1, s=50, 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(lon_distusavaf2,lat_distusavaf2)
csf2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colsavaf,                 alpha=0.8, s=50, zorder=4, label = "Open Woodland / Forest")


#fore
m1_sub_lons, m1_sub_lats = m1(lon_distufore2,lat_distufore2)
cfs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colfore,                 alpha=0.8, s=50, zorder=4, label = "Forest / Open Woodland")

# Shrink current axis by 20%
#box = ax.get_position()
#ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
legend1 = plt.legend((cbs, csb, csf, cfs), ('Grassland / Open Woodland', 'Open Woodland / Grassland', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1., 0.4), markerscale = 2, ncol = 1, fontsize = 14, fancybox = True, frameon=True, title = "Fire Disturbed Tree Cover States")

legend2 = plt.legend((cbs2, csb2, csf2, cfs2), ('Grassland / Open Woodland', 'Open Woodland / Grassland', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1., 0.6), markerscale = 2, ncol = 1, fontsize = 14, fancybox = True, frameon=True, title = "Equivalent Tree Cover States")
legend2.get_title().set_fontsize('12') #legend 'Title' fontsize
legend1.get_title().set_fontsize('12') #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/Nea_East_Possible_Alternative_001.pdf', dpi=600)


plt.show()


# In[14]:

# plot it [white plot] # states where all variables are the same but without repetition with fire stuff
fig = plt.figure(figsize=(12,12), facecolor='white')
# use generative function for legend, requires a figure, 4 data_arrays, 4 lists for legend entries, optionally list of 4 cmaps
#el_bare = [el for el in np.unique(c_distubare2) if np.ma.is_masked(el)==False]
#el_savab = [el for el in np.unique(c_distusavab2) if np.ma.is_masked(el)==False]
#el_savaf = [el for el in np.unique(c_distusavaf2) if np.ma.is_masked(el)==False]
#el_fore = [el for el in np.unique(c_distufore2) if np.ma.is_masked(el)==False]
#patches = generate_legend(fig,c_distubare2, c_distusavab2, c_distufore2, c_distusavaf2, el_bare, el_savab, el_fore, el_savaf, cmaps=[cmapbare_r, cmapsavab_r, cmapfore_r, cmapsava_r])
# create ax and set title
ax = fig.add_subplot(111)
ax.set_title("Same Variables Tree Cover States over Northern East Eurasia", fontsize=18)
#m = Basemap(projection='ortho',lon_0=-125,lat_0=60,resolution='l') #West
m = Basemap(projection='ortho',lon_0=100.,lat_0=60,resolution='l') #East
width = m.urcrnrx - m.llcrnrx
height = m.urcrnry - m.llcrnry
# let's zoom in the plot, cause the original basemap is too zoomed out
coef = 0.45 # zooming factor (percentage of original dimensions)
width = width*coef
height = height*coef
m1 = Basemap(projection='ortho',lon_0=100.,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
par = m1.drawparallels(np.arange(-80.,81.,20.),labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,181.,20.),labels=[False,False,False,True], color='black', fontsize=14)
m1.drawmapboundary(fill_color='black')
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')

#bare
m1_sub_lons, m1_sub_lats = m1(lon_distubare2,lat_distubare2) # convert lon lat to basemap coordinates
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colbare,                 alpha=0.8, s=50, 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(lon_distusavab2,lat_distusavab2)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colsavab,                 alpha=1, s=50, 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(lon_distusavaf2,lat_distusavaf2)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colsavaf,                 alpha=0.8, s=50, zorder=4, label = "Open Woodland / Forest")


#fore
m1_sub_lons, m1_sub_lats = m1(lon_distufore2,lat_distufore2)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colfore,                 alpha=0.8, s=50, zorder=4, label = "Forest / Open Woodland")

legend = plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.03), markerscale = 2, ncol = 2, fontsize = 14, fancybox = True, frameon=True)
#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Meeting_Victor/Alternative_States/New_Classification_with_SoilTexture_and_Disturbances/Nea_East_Equivalent_States_001.pdf', dpi=600)

plt.show()

# plot it [white plot] # states where all variables are the same but without repetition with fire stuff
fig = plt.figure(figsize=(12,12), facecolor='white')
# use generative function for legend, requires a figure, 4 data_arrays, 4 lists for legend entries, optionally list of 4 cmaps
#el_bare = [el for el in np.unique(c_distubare2) if np.ma.is_masked(el)==False]
#el_savab = [el for el in np.unique(c_distusavab2) if np.ma.is_masked(el)==False]
#el_savaf = [el for el in np.unique(c_distusavaf2) if np.ma.is_masked(el)==False]
#el_fore = [el for el in np.unique(c_distufore2) if np.ma.is_masked(el)==False]
#patches = generate_legend(fig,c_distubare2, c_distusavab2, c_distufore2, c_distusavaf2, el_bare, el_savab, el_fore, el_savaf, cmaps=[cmapbare_r, cmapsavab_r, cmapfore_r, cmapsava_r])
# create ax and set title
ax = fig.add_subplot(111)
ax.set_title("Same Variables Tree Cover States over Northern East Eurasia", fontsize=18)
#m = Basemap(projection='ortho',lon_0=-125,lat_0=60,resolution='l') #West
m = Basemap(projection='ortho',lon_0=100.,lat_0=60,resolution='l') #East
width = m.urcrnrx - m.llcrnrx
height = m.urcrnry - m.llcrnry
# let's zoom in the plot, cause the original basemap is too zoomed out
coef = 0.45 # zooming factor (percentage of original dimensions)
width = width*coef
height = height*coef
m1 = Basemap(projection='ortho',lon_0=100.,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
par = m1.drawparallels(np.arange(-80.,81.,20.),labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,181.,20.),labels=[False,False,False,True], color='black', fontsize=14)
m1.drawmapboundary(fill_color='black')
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')

#bare
m1_sub_lons, m1_sub_lats = m1(lon_distubare2,lat_distubare2) # convert lon lat to basemap coordinates
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=50, 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(lon_distusavab2,lat_distusavab2)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=1, s=50, 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(lon_distusavaf2,lat_distusavaf2)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=50, zorder=4, label = "Open Woodland / Forest")


#fore
m1_sub_lons, m1_sub_lats = m1(lon_distufore2,lat_distufore2)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=50, zorder=4, label = "Forest / Open Woodland")

legend = plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.03), markerscale = 2, ncol = 2, fontsize = 14, fancybox = True, frameon=True)



plt.show()


# In[36]:

# possible states, but with different projection
fig = plt.figure(figsize=(12,12), facecolor='white')
# use generative function for legend, requires a figure, 4 data_arrays, 4 lists for legend entries, optionally list of 4 cmaps
#el_bare = [el for el in np.unique(c_distubare2) if np.ma.is_masked(el)==False]
#el_savab = [el for el in np.unique(c_distusavab2) if np.ma.is_masked(el)==False]
#el_savaf = [el for el in np.unique(c_distusavaf2) if np.ma.is_masked(el)==False]
#el_fore = [el for el in np.unique(c_distufore2) if np.ma.is_masked(el)==False]
#patches = generate_legend(fig,c_distubare2, c_distusavab2, c_distufore2, c_distusavaf2, el_bare, el_savab, el_fore, el_savaf, cmaps=[cmapbare_r, cmapsavab_r, cmapfore_r, cmapsava_r])
# create ax and set title
ax = fig.add_subplot(111)
ax.set_title("Same Variables Tree Cover States over Northern East Eurasia", fontsize=18)
#m = Basemap(projection='ortho',lon_0=-125,lat_0=60,resolution='l') #West
m = Basemap(projection='ortho',lon_0=100.,lat_0=60,resolution='l') #East
width = m.urcrnrx - m.llcrnrx
height = m.urcrnry - m.llcrnry
# let's zoom in the plot, cause the original basemap is too zoomed out
coef = 0.45 # zooming factor (percentage of original dimensions)
width = width*coef
height = height*coef
m1 = Basemap(projection='ortho',lon_0=100.,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
par = m1.drawparallels(np.arange(-80.,81.,20.),labels=[True,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,181.,20.),labels=[False,False,False,True], color='black', fontsize=14)
m1.drawmapboundary(fill_color='black')
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')


#bare1
m1_sub_lons, m1_sub_lats = m1(lon_bare_1_sb,lat_bare_1_sb)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=20, zorder=3, label = 'Grassland / Open Woodland')
#bare2
m1_sub_lons, m1_sub_lats = m1(lon_bare_2_bs,lat_bare_2_bs)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colbare,                 alpha=0.8, s=20, zorder=3)
#savab1
m1_sub_lons, m1_sub_lats = m1(lon_sava_1_sb,lat_sava_1_sb)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=0.8, s=20, zorder=4, label = 'Open Woodland / Grassland')
#savab2
m1_sub_lons, m1_sub_lats = m1(lon_sava_2_bs,lat_sava_2_bs)
csb = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsavab,                 alpha=0.8, s=20, zorder=4)
#savaf1
m1_sub_lons, m1_sub_lats = m1(lon_sava_1_sf,lat_sava_1_sf)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=20, zorder=4, label = 'Open Woodland / Forest')
#savaf2
m1_sub_lons, m1_sub_lats = m1(lon_sava_2_fs,lat_sava_2_fs)
csf = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colsava,                 alpha=0.8, s=20, zorder=4)
#fore1
m1_sub_lons, m1_sub_lats = m1(lon_fore_1_sf,lat_fore_1_sf)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=20, zorder=4, label = 'Forest / Open Woodland')
#fore2
m1_sub_lons, m1_sub_lats = m1(lon_fore_2_fs,lat_fore_2_fs)
cfs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=colfore,                 alpha=0.8, s=20, zorder=4)


# not fire disturbed
#bare
m1_sub_lons, m1_sub_lats = m1(lon_distubare2,lat_distubare2) # convert lon lat to basemap coordinates
cbs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colbare,                 alpha=0.8, s=20, 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(lon_distusavab2,lat_distusavab2)
csb2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colsavab,                 alpha=1, s=20, 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(lon_distusavaf2,lat_distusavaf2)
csf2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colsavaf,                 alpha=0.8, s=20, zorder=4, label = "Open Woodland / Forest")


#fore
m1_sub_lons, m1_sub_lats = m1(lon_distufore2,lat_distufore2)
cfs2 = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=alt_colfore,                 alpha=0.8, s=20, zorder=4, label = "Forest / Open Woodland")

# Shrink current axis by 20%
#box = ax.get_position()
#ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
legend1 = plt.legend((cbs, csb, csf, cfs), ('Grassland / Open Woodland', 'Open Woodland / Grassland', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1., 0.4), markerscale = 2, ncol = 1, fontsize = 14, fancybox = True, frameon=True,                      title = "Fire Disturbed Tree Cover States")

legend2 = plt.legend((cbs2, csb2, csf2, cfs2), ('Grassland / Open Woodland', 'Open Woodland / Grassland', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1., 0.6), markerscale = 2, ncol = 1, fontsize = 14, fancybox = True, frameon=True,                      title = "Equivalent Tree Cover States")
legend2.get_title().set_fontsize('12') #legend 'Title' fontsize
legend1.get_title().set_fontsize('12') #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/Nea_East_Possible_Alternative_001.pdf', dpi=600)


plt.show()


# In[60]:

# possible states, but with different projection
fig = plt.figure(figsize=(15,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
par = m1.drawparallels([30.,45.,55.,65.,75.,90],labels=[False,False,False,False], color='black', fontsize=14)
mer = m1.drawmeridians(np.arange(-180.,179.,45.),labels=[False,False,True,True], color='black', fontsize=14)
#list of 4 values (default [0,0,0,0]) that control
#whether parallels are labelled where they intersect
#the left, right, top or bottom of the plot. For
#example labels=[1,0,0,1] will cause parallels
#to be labelled where they intersect the left and
#and bottom of the plot, but not the right and top.
#labelstyle       if set to "+/-", north and south latitudes are
#labelled with "+" and "-", otherwise they are
#labelled with "N" and "S".
m1.drawmapboundary(fill_color='#02d8e9')#black
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
altst_plot(m1=m1, lonlat = ea_e_dist_lonlat, dotsize = 20, alpha=0.8, marker = 'h')

# nondist
st_plot(m1=m1, lonlat = ea_e_alt1_lonlat, dotsize = 20, alpha=0.8, marker = 'h')

st_plot(m1=m1, lonlat = ea_e_alt2_lonlat, dotsize = 20, alpha=0.8, marker = 'h')


# Shrink current axis by 20%
#box = ax.get_position()
#ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
legend1 = plt.legend((cbs, csb, csf, cfs), ('Grassland / Open Woodland', 'Open Woodland / Grassland', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1.05, 0.4), markerscale = 2, ncol = 1, fontsize = 14, fancybox = True, frameon=True,                      title = "Fire Disturbed Tree Cover States")

legend2 = plt.legend((cbs2, csb2, csf2, cfs2), ('Grassland / Open Woodland', 'Open Woodland / Grassland', 'Open Woodland / Forest', 'Forest / Open Woodland') ,                     loc='center left', bbox_to_anchor=(1.05, 0.6), markerscale = 2, ncol = 1, fontsize = 14, fancybox = True, frameon=True,                      title = "Equivalent Tree Cover States")
legend2.get_title().set_fontsize('12') #legend 'Title' fontsize
legend1.get_title().set_fontsize('12') #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/Nea_East_Possible_Alternative_001.pdf', dpi=600)


plt.show()

