
# coding: utf-8

# In[5]:

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import seaborn as sns
import datetime
import rpy2.robjects as robjects
import rpy2.robjects.packages as rpackages # import rpy2's package module
from rpy2.robjects.packages import importr # import R's utility package
import netCDF4
import sys
import math 
import numpy.ma as ma
from matplotlib.colors import ListedColormap

# sys.path.append("/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/") # Edited .zshrc profile with PYTHONPATH 
from Functions import plot_north_polar
from pandas import DataFrame 
from collections import Counter


r = robjects.r
# Some R packages
mgcv = importr('mgcv')
importr("data.table")
utils = importr('utils')
# select a mirror for R packages
utils.chooseCRANmirror(ind=1) # select the first mirror in the list
from rpy2.robjects import FloatVector
import rpy2.robjects.numpy2ri # For the conversion of numpy objects into R 
rpy2.robjects.numpy2ri.activate()
silvermantest = importr('silvermantest')


sns.set_style(style="white")






#final arrays with lon, lat of all alternative states 
commonlons_bor = np.loadtxt('commonlons_bor.txt')
commonlats_bor = np.loadtxt('commonlats_bor.txt')
commonlai_bor = np.loadtxt('commonlai_bor.txt')
commonndvi_bor = np.loadtxt('commonndvi_bor.txt')
commontcf_bor = np.loadtxt('commontcf_bor.txt')
commonstate_bor = np.loadtxt('commonstate_bor.txt')
orderednums = np.loadtxt('orderednums.txt')


# let's divide TCF into the three states
tcf_less = [elem for elem in commontcf_bor if elem<20]
tcf_open = [elem for elem in commontcf_bor if (elem>=20 and elem<45)]
tcf_fore = [elem for elem in commontcf_bor if elem>=45]




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

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


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)




# In[6]:

# execute this only if you want to re compute the variables

path = '/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Boreal_New_Classification/'

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


# concatenate lonlat arrays of fire disturbed and equivalent tree cover arrays
# create same shape arrays filled with increasing numbers
altlist = [na_e_dist_lonlat, na_w_dist_lonlat, ea_e_dist_lonlat, ea_w_dist_lonlat]
altlon = np.concatenate([np.concatenate([ma.compressed(elem[i*2]) for i in range(4)]) for elem in altlist])
altlat = np.concatenate([np.concatenate([ma.compressed(elem[i*2+1]) for i in range(4)]) for elem in altlist])
altnum = np.concatenate([np.concatenate([np.ones(ma.compressed(elem[i*2+1]).shape)*i for i in range(4)]) for elem in altlist])

fdlist = [na_e_alt1_lonlat, na_w_alt1_lonlat, ea_e_alt1_lonlat, ea_w_alt1_lonlat,            na_e_alt2_lonlat, na_w_alt2_lonlat, ea_e_alt2_lonlat, ea_w_alt2_lonlat] 
fdlon = np.concatenate([np.concatenate([ma.compressed(elem[i*2]) for i in range(4)]) for elem in fdlist])
fdlat = np.concatenate([np.concatenate([ma.compressed(elem[i*2+1]) for i in range(4)]) for elem in fdlist])
fdnum = np.concatenate([np.concatenate([np.ones(ma.compressed(elem[i*2+1]).shape)*(i+4) for i in range(4)]) for elem in fdlist])

# concatenate into single arrays to have everything together
altfdlon = np.concatenate([altlon,fdlon])
altfdlat = np.concatenate([altlat,fdlat])
altfdnum = np.concatenate([altnum,fdnum])

#create cmap for alternative states
colseq = [colbare,colsavab,colsava,colfore,alt_colbare,alt_colsavab,alt_colsavaf,alt_colfore]
current_palette = sns.color_palette(colseq)
altfdmap = ListedColormap(current_palette)

# sort a lons of the altfd array so that they match the commonlons_bor, and same for the numerical array necessary for colouring
objmap = {altfdlon[i]: altfdnum[i] for i in range(len(altfdlon))}
orderednums = [objmap[id] for id in commonlons_bor]
#np.savetxt('orderednums.txt', orderednums)


# In[7]:

# let's check bimodality of every combination (open treeless, open forest), for every region->next put the regions together


# In[8]:

# creates the variables for comparisons
import pandas as pd
zone = '_ea_east'
path_to_pickles = '/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/'
bare_classified_ea_e = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')
savanna_classified_ea_e = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')
forest_classified_ea_e = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')

treeless_open_common_ea_e = np.array([elem for elem in np.array(bare_classified_ea_e) if np.any(elem[0]==np.array(savanna_classified_ea_e)[:,0])])
open_treeless_common_ea_e = np.array([elem for elem in np.array(savanna_classified_ea_e) if np.any(elem[0]==np.array(bare_classified_ea_e)[:,0])])

op_less_common_ea_e = np.concatenate([treeless_open_common_ea_e, open_treeless_common_ea_e])
op_less_lat_ea_e, op_less_lon_ea_e, op_less_state_ea_e = op_less_common_ea_e[:, -1], op_less_common_ea_e[:, -2], op_less_common_ea_e[:, 0]
common_op_less_state_ea_e = np.array([op_less_state_ea_e[(op_less_lon_ea_e==commonlons_bor[ind]) & (op_less_lat_ea_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_ea_e[(op_less_lon_ea_e==commonlons_bor[ind]) & (op_less_lat_ea_e==commonlats_bor[ind])])>0)])
common_op_less_state_lon_ea_e = np.array([op_less_lon_ea_e[(op_less_lon_ea_e==commonlons_bor[ind]) & (op_less_lat_ea_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_ea_e[(op_less_lon_ea_e==commonlons_bor[ind]) & (op_less_lat_ea_e==commonlats_bor[ind])])>0)])
common_op_less_state_lat_ea_e = np.array([op_less_lat_ea_e[(op_less_lon_ea_e==commonlons_bor[ind]) & (op_less_lat_ea_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_ea_e[(op_less_lon_ea_e==commonlons_bor[ind]) & (op_less_lat_ea_e==commonlats_bor[ind])])>0)])
bare_open_tcf_comparison_ea_e = np.array([commontcf_bor[((commonlons_bor==common_op_less_state_lon_ea_e[ind]) & (commonlats_bor==common_op_less_state_lat_ea_e[ind]))][0] for ind in range(len(common_op_less_state_lat_ea_e))])


forest_open_common_ea_e = np.array([elem for elem in np.array(forest_classified_ea_e) if np.any(elem[0]==np.array(savanna_classified_ea_e)[:,0])])
open_forest_common_ea_e = np.array([elem for elem in np.array(savanna_classified_ea_e) if np.any(elem[0]==np.array(forest_classified_ea_e)[:,0])])

op_fore_common_ea_e = np.concatenate([forest_open_common_ea_e, open_forest_common_ea_e])
op_fore_lat_ea_e, op_fore_lon_ea_e, op_fore_state_ea_e = op_fore_common_ea_e[:, -1], op_fore_common_ea_e[:, -2], op_fore_common_ea_e[:, 0]
common_op_fore_state_ea_e = np.array([op_fore_state_ea_e[(op_fore_lon_ea_e==commonlons_bor[ind]) & (op_fore_lat_ea_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_ea_e[(op_fore_lon_ea_e==commonlons_bor[ind]) & (op_fore_lat_ea_e==commonlats_bor[ind])])>0)])
common_op_fore_state_lon_ea_e = np.array([op_fore_lon_ea_e[(op_fore_lon_ea_e==commonlons_bor[ind]) & (op_fore_lat_ea_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_ea_e[(op_fore_lon_ea_e==commonlons_bor[ind]) & (op_fore_lat_ea_e==commonlats_bor[ind])])>0)])
common_op_fore_state_lat_ea_e = np.array([op_fore_lat_ea_e[(op_fore_lon_ea_e==commonlons_bor[ind]) & (op_fore_lat_ea_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_ea_e[(op_fore_lon_ea_e==commonlons_bor[ind]) & (op_fore_lat_ea_e==commonlats_bor[ind])])>0)])
forest_open_tcf_comparison_ea_e = np.array([commontcf_bor[((commonlons_bor==common_op_fore_state_lon_ea_e[ind]) & (commonlats_bor==common_op_fore_state_lat_ea_e[ind]))][0] for ind in range(len(common_op_fore_state_lat_ea_e))])

#r.assign('tcf_bor_op_less_ea_e', np.array(bare_open_tcf_comparison_ea_e))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_ea_e),1, adjust=TRUE)');
#print r('stest')

#r.assign('tcf_bor_op_fore_ea_e', np.array(forest_open_tcf_comparison_ea_e))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_ea_e),1, adjust=TRUE)');
#print r('stest')

#EA_W
zone = '_ea_west_russia'
bare_classified_ea_w = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')
savanna_classified_ea_w = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')
forest_classified_ea_w = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')

treeless_open_common_ea_w = np.array([elem for elem in np.array(bare_classified_ea_w) if np.any(elem[0]==np.array(savanna_classified_ea_w)[:,0])])
open_treeless_common_ea_w = np.array([elem for elem in np.array(savanna_classified_ea_w) if np.any(elem[0]==np.array(bare_classified_ea_w)[:,0])])

op_less_common_ea_w = np.concatenate([treeless_open_common_ea_w, open_treeless_common_ea_w])
op_less_lat_ea_w, op_less_lon_ea_w, op_less_state_ea_w = op_less_common_ea_w[:, -1], op_less_common_ea_w[:, -2], op_less_common_ea_w[:, 0]
common_op_less_state_ea_w = np.array([op_less_state_ea_w[(op_less_lon_ea_w==commonlons_bor[ind]) & (op_less_lat_ea_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_ea_w[(op_less_lon_ea_w==commonlons_bor[ind]) & (op_less_lat_ea_w==commonlats_bor[ind])])>0)])
common_op_less_state_lon_ea_w = np.array([op_less_lon_ea_w[(op_less_lon_ea_w==commonlons_bor[ind]) & (op_less_lat_ea_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_ea_w[(op_less_lon_ea_w==commonlons_bor[ind]) & (op_less_lat_ea_w==commonlats_bor[ind])])>0)])
common_op_less_state_lat_ea_w = np.array([op_less_lat_ea_w[(op_less_lon_ea_w==commonlons_bor[ind]) & (op_less_lat_ea_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_ea_w[(op_less_lon_ea_w==commonlons_bor[ind]) & (op_less_lat_ea_w==commonlats_bor[ind])])>0)])
bare_open_tcf_comparison_ea_w = np.array([commontcf_bor[((commonlons_bor==common_op_less_state_lon_ea_w[ind]) & (commonlats_bor==common_op_less_state_lat_ea_w[ind]))][0] for ind in range(len(common_op_less_state_lat_ea_w))])


forest_open_common_ea_w = np.array([elem for elem in np.array(forest_classified_ea_w) if np.any(elem[0]==np.array(savanna_classified_ea_w)[:,0])])
open_forest_common_ea_w = np.array([elem for elem in np.array(savanna_classified_ea_w) if np.any(elem[0]==np.array(forest_classified_ea_w)[:,0])])

op_fore_common_ea_w = np.concatenate([forest_open_common_ea_w, open_forest_common_ea_w])
op_fore_lat_ea_w, op_fore_lon_ea_w, op_fore_state_ea_w = op_fore_common_ea_w[:, -1], op_fore_common_ea_w[:, -2], op_fore_common_ea_w[:, 0]
common_op_fore_state_ea_w = np.array([op_fore_state_ea_w[(op_fore_lon_ea_w==commonlons_bor[ind]) & (op_fore_lat_ea_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_ea_w[(op_fore_lon_ea_w==commonlons_bor[ind]) & (op_fore_lat_ea_w==commonlats_bor[ind])])>0)])
common_op_fore_state_lon_ea_w = np.array([op_fore_lon_ea_w[(op_fore_lon_ea_w==commonlons_bor[ind]) & (op_fore_lat_ea_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_ea_w[(op_fore_lon_ea_w==commonlons_bor[ind]) & (op_fore_lat_ea_w==commonlats_bor[ind])])>0)])
common_op_fore_state_lat_ea_w = np.array([op_fore_lat_ea_w[(op_fore_lon_ea_w==commonlons_bor[ind]) & (op_fore_lat_ea_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_ea_w[(op_fore_lon_ea_w==commonlons_bor[ind]) & (op_fore_lat_ea_w==commonlats_bor[ind])])>0)])
forest_open_tcf_comparison_ea_w = np.array([commontcf_bor[((commonlons_bor==common_op_fore_state_lon_ea_w[ind]) & (commonlats_bor==common_op_fore_state_lat_ea_w[ind]))][0] for ind in range(len(common_op_fore_state_lat_ea_w))])

#r.assign('tcf_bor_op_less_ea_w', np.array(bare_open_tcf_comparison_ea_w))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_ea_w),1, adjust=TRUE)');
#print r('stest')

#r.assign('tcf_bor_op_fore_ea_w', np.array(forest_open_tcf_comparison_ea_w))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_ea_w),1, adjust=TRUE)');
#print r('stest')


# NA_E
zone = '_na_east'
bare_classified_na_e = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')
savanna_classified_na_e = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')
forest_classified_na_e = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')

treeless_open_common_na_e = np.array([elem for elem in np.array(bare_classified_na_e) if np.any(elem[0]==np.array(savanna_classified_na_e)[:,0])])
open_treeless_common_na_e = np.array([elem for elem in np.array(savanna_classified_na_e) if np.any(elem[0]==np.array(bare_classified_na_e)[:,0])])

op_less_common_na_e = np.concatenate([treeless_open_common_na_e, open_treeless_common_na_e])
op_less_lat_na_e, op_less_lon_na_e, op_less_state_na_e = op_less_common_na_e[:, -1], op_less_common_na_e[:, -2], op_less_common_na_e[:, 0]
common_op_less_state_na_e = np.array([op_less_state_na_e[(op_less_lon_na_e==commonlons_bor[ind]) & (op_less_lat_na_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_na_e[(op_less_lon_na_e==commonlons_bor[ind]) & (op_less_lat_na_e==commonlats_bor[ind])])>0)])
common_op_less_state_lon_na_e = np.array([op_less_lon_na_e[(op_less_lon_na_e==commonlons_bor[ind]) & (op_less_lat_na_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_na_e[(op_less_lon_na_e==commonlons_bor[ind]) & (op_less_lat_na_e==commonlats_bor[ind])])>0)])
common_op_less_state_lat_na_e = np.array([op_less_lat_na_e[(op_less_lon_na_e==commonlons_bor[ind]) & (op_less_lat_na_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_na_e[(op_less_lon_na_e==commonlons_bor[ind]) & (op_less_lat_na_e==commonlats_bor[ind])])>0)])
bare_open_tcf_comparison_na_e = np.array([commontcf_bor[((commonlons_bor==common_op_less_state_lon_na_e[ind]) & (commonlats_bor==common_op_less_state_lat_na_e[ind]))][0] for ind in range(len(common_op_less_state_lat_na_e))])


forest_open_common_na_e = np.array([elem for elem in np.array(forest_classified_na_e) if np.any(elem[0]==np.array(savanna_classified_na_e)[:,0])])
open_forest_common_na_e = np.array([elem for elem in np.array(savanna_classified_na_e) if np.any(elem[0]==np.array(forest_classified_na_e)[:,0])])

op_fore_common_na_e = np.concatenate([forest_open_common_na_e, open_forest_common_na_e])
op_fore_lat_na_e, op_fore_lon_na_e, op_fore_state_na_e = op_fore_common_na_e[:, -1], op_fore_common_na_e[:, -2], op_fore_common_na_e[:, 0]
common_op_fore_state_na_e = np.array([op_fore_state_na_e[(op_fore_lon_na_e==commonlons_bor[ind]) & (op_fore_lat_na_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_na_e[(op_fore_lon_na_e==commonlons_bor[ind]) & (op_fore_lat_na_e==commonlats_bor[ind])])>0)])
common_op_fore_state_lon_na_e = np.array([op_fore_lon_na_e[(op_fore_lon_na_e==commonlons_bor[ind]) & (op_fore_lat_na_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_na_e[(op_fore_lon_na_e==commonlons_bor[ind]) & (op_fore_lat_na_e==commonlats_bor[ind])])>0)])
common_op_fore_state_lat_na_e = np.array([op_fore_lat_na_e[(op_fore_lon_na_e==commonlons_bor[ind]) & (op_fore_lat_na_e==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_na_e[(op_fore_lon_na_e==commonlons_bor[ind]) & (op_fore_lat_na_e==commonlats_bor[ind])])>0)])
forest_open_tcf_comparison_na_e = np.array([commontcf_bor[((commonlons_bor==common_op_fore_state_lon_na_e[ind]) & (commonlats_bor==common_op_fore_state_lat_na_e[ind]))][0] for ind in range(len(common_op_fore_state_lat_na_e))])

#r.assign('tcf_bor_op_less_na_e', np.array(bare_open_tcf_comparison_na_e))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_na_e),1, adjust=TRUE)');
#print r('stest')

#r.assign('tcf_bor_op_fore_na_e', np.array(forest_open_tcf_comparison_na_e))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_na_e),1, adjust=TRUE)');
#print r('stest')

#NA_W
zone = '_na_west'
bare_classified_na_w = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')
savanna_classified_na_w = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')
forest_classified_na_w = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')

treeless_open_common_na_w = np.array([elem for elem in np.array(bare_classified_na_w) if np.any(elem[0]==np.array(savanna_classified_na_w)[:,0])])
open_treeless_common_na_w = np.array([elem for elem in np.array(savanna_classified_na_w) if np.any(elem[0]==np.array(bare_classified_na_w)[:,0])])

op_less_common_na_w = np.concatenate([treeless_open_common_na_w, open_treeless_common_na_w])
op_less_lat_na_w, op_less_lon_na_w, op_less_state_na_w = op_less_common_na_w[:, -1], op_less_common_na_w[:, -2], op_less_common_na_w[:, 0]
common_op_less_state_na_w = np.array([op_less_state_na_w[(op_less_lon_na_w==commonlons_bor[ind]) & (op_less_lat_na_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_na_w[(op_less_lon_na_w==commonlons_bor[ind]) & (op_less_lat_na_w==commonlats_bor[ind])])>0)])
common_op_less_state_lon_na_w = np.array([op_less_lon_na_w[(op_less_lon_na_w==commonlons_bor[ind]) & (op_less_lat_na_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_na_w[(op_less_lon_na_w==commonlons_bor[ind]) & (op_less_lat_na_w==commonlats_bor[ind])])>0)])
common_op_less_state_lat_na_w = np.array([op_less_lat_na_w[(op_less_lon_na_w==commonlons_bor[ind]) & (op_less_lat_na_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_less_state_na_w[(op_less_lon_na_w==commonlons_bor[ind]) & (op_less_lat_na_w==commonlats_bor[ind])])>0)])
bare_open_tcf_comparison_na_w = np.array([commontcf_bor[((commonlons_bor==common_op_less_state_lon_na_w[ind]) & (commonlats_bor==common_op_less_state_lat_na_w[ind]))][0] for ind in range(len(common_op_less_state_lat_na_w))])


forest_open_common_na_w = np.array([elem for elem in np.array(forest_classified_na_w) if np.any(elem[0]==np.array(savanna_classified_na_w)[:,0])])
open_forest_common_na_w = np.array([elem for elem in np.array(savanna_classified_na_w) if np.any(elem[0]==np.array(forest_classified_na_w)[:,0])])

op_fore_common_na_w = np.concatenate([forest_open_common_na_w, open_forest_common_na_w])
op_fore_lat_na_w, op_fore_lon_na_w, op_fore_state_na_w = op_fore_common_na_w[:, -1], op_fore_common_na_w[:, -2], op_fore_common_na_w[:, 0]
common_op_fore_state_na_w = np.array([op_fore_state_na_w[(op_fore_lon_na_w==commonlons_bor[ind]) & (op_fore_lat_na_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_na_w[(op_fore_lon_na_w==commonlons_bor[ind]) & (op_fore_lat_na_w==commonlats_bor[ind])])>0)])
common_op_fore_state_lon_na_w = np.array([op_fore_lon_na_w[(op_fore_lon_na_w==commonlons_bor[ind]) & (op_fore_lat_na_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_na_w[(op_fore_lon_na_w==commonlons_bor[ind]) & (op_fore_lat_na_w==commonlats_bor[ind])])>0)])
common_op_fore_state_lat_na_w = np.array([op_fore_lat_na_w[(op_fore_lon_na_w==commonlons_bor[ind]) & (op_fore_lat_na_w==commonlats_bor[ind])][0] for ind in range(len(commonlats_bor)) if (len(op_fore_state_na_w[(op_fore_lon_na_w==commonlons_bor[ind]) & (op_fore_lat_na_w==commonlats_bor[ind])])>0)])
forest_open_tcf_comparison_na_w = np.array([commontcf_bor[((commonlons_bor==common_op_fore_state_lon_na_w[ind]) & (commonlats_bor==common_op_fore_state_lat_na_w[ind]))][0] for ind in range(len(common_op_fore_state_lat_na_w))])

# not enough points to compute
#r.assign('tcf_bor_op_less_na_w', np.array(bare_open_tcf_comparison_na_w))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_na_w),1, adjust=TRUE)');
#print r('stest')

#r.assign('tcf_bor_op_fore_na_w', np.array(forest_open_tcf_comparison_na_w))
#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_na_w),1, adjust=TRUE)');
#print r('stest')



# lon lat for boreal open treeless and open forest

forest_open_bor_lon = np.concatenate([common_op_fore_state_lon_ea_e, common_op_fore_state_lon_ea_w, common_op_fore_state_lon_na_e, common_op_fore_state_lon_na_w])
forest_open_bor_lat = np.concatenate([common_op_fore_state_lat_ea_e, common_op_fore_state_lat_ea_w, common_op_fore_state_lat_na_e, common_op_fore_state_lat_na_w])

open_treeless_bor_lon = np.concatenate([common_op_less_state_lon_ea_e, common_op_less_state_lon_ea_w, common_op_less_state_lon_na_e, common_op_less_state_lon_na_w])
open_treeless_bor_lat = np.concatenate([common_op_less_state_lat_ea_e, common_op_less_state_lat_ea_w, common_op_less_state_lat_na_e, common_op_less_state_lat_na_w])


#entire boreal area
open_treeless_bor = np.concatenate([bare_open_tcf_comparison_ea_e, bare_open_tcf_comparison_ea_w, bare_open_tcf_comparison_na_e, bare_open_tcf_comparison_na_w])
forest_open_bor = np.concatenate([forest_open_tcf_comparison_ea_e, forest_open_tcf_comparison_ea_w, forest_open_tcf_comparison_na_e, forest_open_tcf_comparison_na_w])

#open_treeless_bor = np.loadtxt("open_treless_bor.txt")
#forest_open_bor = np.loadtxt("forest_open_bor.txt")
# let's divide them into the single states

alt_less = [elem for elem in open_treeless_bor if elem<20]
alt_open_l = [elem for elem in open_treeless_bor if (elem>=20)]
alt_open_f = [elem for elem in forest_open_bor if (elem<45)]
alt_fore = [elem for elem in forest_open_bor if elem>=45]

# Fitting of the histograms with KDE
from scipy import stats
import statsmodels.api as sm
from statsmodels.distributions.mixture_rvs import mixture_rvs
from statsmodels.nonparametric.kernel_density import KDEMultivariate as kdem
import peakutils
kde = sm.nonparametric.KDEUnivariate(open_treeless_bor)
kde.fit(bw=4)
kde2 = sm.nonparametric.KDEUnivariate(forest_open_bor)
kde2.fit(bw=4)

indexes = peakutils.indexes(kde.density, thres=0.02)#, thres=0.02/max(cb), min_dist=100)
indexes2 = peakutils.indexes(kde2.density, thres=0.02)#, thres=0.02/max(cb), min_dist=100)
peaks_x = [kde.support[i] for i in indexes] 
peaks_y = [kde.density[i] for i in indexes]
peaks_x2 = [kde2.support[i] for i in indexes2] 
peaks_y2 = [kde2.density[i] for i in indexes2]



min_peak_dist = np.min([peaks_x[1]-peaks_x[0], peaks_x2[1]-peaks_x2[0]])
print "minimum peaks distance", min_peak_dist
print "all peaks distances", [peaks_x[1]-peaks_x[0], peaks_x2[1]-peaks_x2[0]]



#Variability

nc_tcf_years = netCDF4.Dataset('/Users/Beniamino/Documents/Datasets/MODIS/MODIS_Aall_Variability/Boreal/Bor_180_MOD44B_V5__VCFFplusVCF__2001-2010__0.5deg__UHAM-ICDC.nc') 

tcf_years = ma.getdata(nc_tcf_years.variables['vcff'][:])
tcf_y_lons, tcf_y_lats = nc_tcf_years.variables['lon'][:], nc_tcf_years.variables['lat'][:]

tcf_y_bor = []
for elem in tcf_years:
    tcf_y_bor.append(np.array([elem[tcf_y_lats==commonlats_bor[i],tcf_y_lons==commonlons_bor[i]][0] for i in range(len(commonlons_bor))]))
tcf_y_arr = np.array(tcf_y_bor)


max_array, min_array = [], []
for i in range(tcf_y_arr.shape[1]):
    max_array.append(np.max(tcf_y_arr[:,i]) if (np.max(tcf_y_arr[:,i])!=255) else np.unique(tcf_y_arr[:,i])[-2])
    min_array.append(np.min(tcf_y_arr[:,i]))
    
    
vari_array = (np.array(max_array)-np.array(min_array))
np.mean(vari_array)

std_var_arr = np.array([np.std(tcf_y_arr[:,i]) for i in range(tcf_y_arr.shape[1])])
print "Number of gridcells with variability >"+str(min_peak_dist), len([elem for elem in std_var_arr if elem>min_peak_dist])
import statistics as st
print "Mean standard deviation", np.mean(std_var_arr), "\nMedian standard deviation", st.median(std_var_arr), "\nStandard deviation of standard deviation", st.stdev(std_var_arr), # is the variability smaller than the distance between the modes of the distributions? Let's find the modes


# In[9]:

# Bimodality tests using silverman's test for the Open woodland + Treeless states 
r.assign('open_treeless_bor', np.array(open_treeless_bor))
r('stest <- silverman.test(as.numeric(open_treeless_bor),1, adjust=TRUE)');
print r('stest')

# Bimodality tests using silverman's test for the Open woodland + Forest states 
r.assign('forest_open_bor', np.array(forest_open_bor))
r('stest <- silverman.test(as.numeric(forest_open_bor),1, adjust=TRUE)');
print r('stest')


# In[10]:

#Silverman's test for trimodality of alternative states 
r.assign('tcf_bor', np.array(commontcf_bor))
r('stest <- silverman.test(as.numeric(tcf_bor),2, adjust=FALSE)');
print r('stest')


# In[11]:

# Some plots


# In[12]:

# Bimodality of alternative states categories: Open woodland + Treeless and Open woodland + Forest
# Histograms and silverman's test for bimodality
fig = plt.figure(figsize=(20,9))
ax = fig.add_subplot(121)
ax.set_title('Alternative Treeless / Open Woodland States', fontsize=16)
hi = plt.hist([alt_less, alt_open_l], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#7570b3', '#d95f02']        , label=["Treeless", "Open woodland"])
sns.despine()
plt.xlim(0,50)
plt.xlabel("Tree Cover Fraction - [%]", fontsize=14)
plt.ylabel("Normalised Frequency", fontsize=14)
plt.legend(fontsize=14,ncol=3, loc='upper center')
ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\leq$'+' 1 \n The P-Value is '+r'$<$'+' 0.001',             xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)

ax = fig.add_subplot(122)
ax.set_title('Alternative Forest / Open Woodland States', fontsize=16)
hi = plt.hist([alt_open_f, alt_fore], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#d95f02', '#1b9e77']        , label=["Open woodland", "Forest"])
sns.despine()
plt.xlim(20,80)
plt.ylim(0,0.12)
plt.xlabel("Tree Cover Fraction - [%]", fontsize=14)
plt.ylabel("Normalised Frequency", fontsize=14)
plt.legend(fontsize=14,ncol=3, loc='upper center')
ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\leq$'+' 1 \n The P-Value is '+r'$<$'+' 0.01 ',             xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)

#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Alternative_States_TCF_Bimodality.pdf', dpi=600, transparent = True)

plt.show()


# In[13]:

# Fitting of the distributions of the two classes of alternative states using KDE, 
# peaks distance and variability comparison

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.set_title("Alternative States TCF KDE", fontsize=16)
ax.hist(open_treeless_bor, bins=10, normed=True, color='#7570b3', alpha=0.7, label = 'Treeless / Open woodland States')
ax.plot(kde.support, kde.density, lw=2, color='black', label = 'KDE estimates');


ax = fig.add_subplot(111)
ax.hist(forest_open_bor, bins=10, normed=True, color='#1b9e77', alpha=0.7, label='Forest / Open woodland States')
ax.plot(kde2.support, kde2.density, lw=2, color='black');
plt.xlabel("Tree Cover Fraction - [%]", fontsize=14)
plt.ylabel("Normalised Frequency", fontsize=14)
sns.despine()
plt.xlim(-10,90)



plt.plot(peaks_x, peaks_y, 'o', color='#d95f02', label = 'KDE peaks')
plt.plot(peaks_x2, peaks_y2, 'o', color='#d95f02')

ax.annotate('Variability of TCF during 2001'+u"\u2013"+'2010: \n Average Standard Deviation is '+str("{:1.2f}".format(np.mean(std_var_arr)))+            '\n Minimum distance between KDE peaks is '+str("{:1.2f}".format(min_peak_dist)),             xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)
plt.legend(fontsize=14)
#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Alternative_States_TCF_Variability.pdf', dpi=600, transparent = True)
plt.show()


# In[ ]:

# Bimodality of alternative states, with Silverman's and parallele to TCF plots
fig = plt.figure(figsize=(15,11))
ax = fig.add_subplot(221)
ax.annotate("Open Woodland / Treeless TCF",  xy=(0.5,1.1), xycoords='axes fraction', ha='center', 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=16)
mer = m1.drawmeridians(np.arange(-135.,179.,45.),labels=[False,False,True,True], color='black', fontsize=16)
m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'
m1.drawcoastlines(linewidth=1.5)
m1.fillcontinents(color='white',lake_color='black')
m1_sub_lons, m1_sub_lats = m1(open_treeless_bor_lon,open_treeless_bor_lat)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=open_treeless_bor, cmap = 'viridis',                 alpha=0.8, s=20, zorder=3)
m1.colorbar(cbs,"right", ticks=np.arange(0,90,10),size="5%", pad="15%")


ax = fig.add_subplot(222)
hi = plt.hist([alt_less, alt_open_l], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#7570b3', '#d95f02']        , label=["Treeless", "Open woodland"])
sns.despine()
plt.xlim(0,50)
plt.xlabel("Tree Cover Fraction - [%]", fontsize=14)
plt.ylabel("Normalised Frequency", fontsize=14)
plt.legend(fontsize=14,ncol=3, loc='upper center')
ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\leq$'+' 1 \n The P-Value is  0.001',             xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)


ax = fig.add_subplot(223)
ax.annotate("Forest / Open Woodland TCF",  xy=(0.5,1.1), xycoords='axes fraction', ha='center', 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=16)
mer = m1.drawmeridians(np.arange(-135.,179.,45.),labels=[False,False,True,True], color='black', fontsize=16)
m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'
m1.drawcoastlines(linewidth=1.5)
m1.fillcontinents(color='white',lake_color='black')
m1_sub_lons, m1_sub_lats = m1(forest_open_bor_lon,forest_open_bor_lat)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=forest_open_bor, cmap = 'viridis',                 alpha=0.8, s=20, zorder=3)
m1.colorbar(cbs,"right", ticks=np.arange(0,90,10),size="5%", pad="15%")


ax = fig.add_subplot(224)
hi = plt.hist([alt_open_f, alt_fore], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#d95f02', '#1b9e77']        , label=["Open woodland", "Forest"])
sns.despine()
plt.xlim(20,80)
plt.ylim(0,0.12)
plt.xlabel("Tree Cover Fraction - [%]", fontsize=14)
plt.ylabel("Normalised Frequency", fontsize=14)
plt.legend(fontsize=14,ncol=3, loc='upper center')
ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\leq$'+' 1 \n The P-Value is 0.01 ',             xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)

plt.show()


# In[ ]:

# TCF plot and histogram for trimodality with silverman's test
fig = plt.figure(figsize=(20,9))
ax = fig.add_subplot(121)
ax.annotate("Alternative States TCF",  xy=(0.5,1.1), xycoords='axes fraction', ha='center', 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=16)
mer = m1.drawmeridians(np.arange(-135.,179.,45.),labels=[False,False,True,True], color='black', fontsize=16)

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


m1_sub_lons, m1_sub_lats = m1(commonlons_bor,commonlats_bor)
cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=commontcf_bor, cmap = 'viridis',                 alpha=0.8, s=20, zorder=3)

m1.colorbar(cbs,"right", ticks=np.arange(0,90,10),size="5%", pad="15%")

ax = fig.add_subplot(122)
ax.set_title('Alternative States TCF Multimodality Test', fontsize=18)
hi = plt.hist([tcf_less, tcf_open, tcf_fore], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#7570b3', '#d95f02', '#1b9e77']        , label=["Treeless", "Open woodland", "Forest"])
sns.despine()
plt.xlim(0,80)
plt.xlabel("Tree Cover Fraction - [%]", fontsize=14)
plt.ylabel("Normalised Frequency", fontsize=14)
plt.legend(fontsize=14,ncol=3, loc='upper center')
ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\leq$'+' 2 \n The P-Value is  0.001001001 ',             xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)

#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Alternative_States_TCF_Multimodality.pdf', dpi=600, transparent = True)

plt.show()


# In[ ]:

# results of silverman's test with different number of modes
stests_ol = [0.002005, 0.3503504, 0.02102102, 0.4614615, 0.4254254] # results pf silvermantests with 1, 2, 3, 4, 5 modes
stests_of = [0.012927, 0.7897898, 0.4744745, 0.4764765, 0.6216216]
x_stests = [1,2,3,4,5]
null_line = np.arange(1,5, 0.01)
x_line = 0.05*np.ones(null_line.shape)
fig = plt.figure()
plt.plot(null_line, x_line, 'r-')
plt.plot(x_stests, stests_ol, '-o', color = '#7570b3', label = 'Treeless / Open woodland')
plt.plot(x_stests, stests_of, '-o', color = '#1b9e77', label = 'Forest / Open woodland')

plt.xlim(0.8,5.2)
plt.ylim(-0.02,1.0)
plt.xlabel('Number of modes in null hypothesis of silvermantest')
plt.ylabel('p-value')
plt.legend()
sns.despine()
plt.show()


# In[ ]:




# In[ ]:

# data load execute to re-create the initial common_variables

#dataforr_ea_west = np.loadtxt("boreal_ea_west_database_lai_ndvi_included.txt")
#data_lons_ea_west = dataforr_ea_west[:,-2]
#data_lats_ea_west = dataforr_ea_west[:,-1]
#data_lai_ea_west = dataforr_ea_west[:,10]
#data_ndvi_ea_west = dataforr_ea_west[:,11]
#data_tcf_ea_west = dataforr_ea_west[:,0]
#data_state_ea_west = dataforr_ea_west[:,-6]



#dataforr_ea_east = np.loadtxt("boreal_ea_east_database_lai_ndvi_included.txt")
#data_lons_ea_east = dataforr_ea_east[:,-2]
#data_lats_ea_east = dataforr_ea_east[:,-1]
#data_lai_ea_east = dataforr_ea_east[:,10]
#data_ndvi_ea_east = dataforr_ea_east[:,11]
#data_tcf_ea_east = dataforr_ea_east[:,0]
#data_state_ea_east = dataforr_ea_east[:,-6]


#dataforr_na_east = np.loadtxt("boreal_na_east_database_lai_ndvi_included.txt")
#data_lons_na_east = dataforr_na_east[:,-2]
#data_lats_na_east = dataforr_na_east[:,-1]
#data_lai_na_east = dataforr_na_east[:,10]
#data_ndvi_na_east = dataforr_na_east[:,11]
#data_tcf_na_east = dataforr_na_east[:,0]
#data_state_na_east = dataforr_na_east[:,-6]#

#dataforr_na_west = np.loadtxt("boreal_na_west_database_lai_ndvi_included.txt")
#data_lons_na_west = dataforr_na_west[:,-2]
#data_lats_na_west = dataforr_na_west[:,-1]
#data_lai_na_west = dataforr_na_west[:,10]
#data_ndvi_na_west = dataforr_na_west[:,11]
#data_tcf_na_west = dataforr_na_west[:,0]
#data_state_na_west = dataforr_na_west[:,-6]#
#
#

## alternative states data#
#

#path = '/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Boreal_New_Classification/'#

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

#ea_w_dist_lons = (ma.concatenate([ea_w_dist_lonlat[0], ea_w_dist_lonlat[2], ea_w_dist_lonlat[4], ea_w_dist_lonlat[6]]))
#ea_w_alt1_lons = (ma.concatenate([ea_w_alt1_lonlat[0], ea_w_alt1_lonlat[2], ea_w_alt1_lonlat[4], ea_w_alt1_lonlat[6]]))
#ea_w_alt2_lons = (ma.concatenate([ea_w_alt2_lonlat[0], ea_w_alt2_lonlat[2], ea_w_alt2_lonlat[4], ea_w_alt2_lonlat[6]]))
#ea_w_dist_lats = (ma.concatenate([ea_w_dist_lonlat[1], ea_w_dist_lonlat[3], ea_w_dist_lonlat[5], ea_w_dist_lonlat[7]]))
#ea_w_alt1_lats = (ma.concatenate([ea_w_alt1_lonlat[1], ea_w_alt1_lonlat[3], ea_w_alt1_lonlat[5], ea_w_alt1_lonlat[7]]))
#ea_w_alt2_lats = (ma.concatenate([ea_w_alt2_lonlat[1], ea_w_alt2_lonlat[3], ea_w_alt2_lonlat[5], ea_w_alt2_lonlat[7]]))
#ea_w_lons = ma.concatenate([ea_w_dist_lons, ea_w_alt1_lons, ea_w_alt2_lons])
#ea_w_lats = ma.concatenate([ea_w_dist_lats, ea_w_alt1_lats, ea_w_alt2_lats])#

#ea_e_dist_lons = (ma.concatenate([ea_e_dist_lonlat[0], ea_e_dist_lonlat[2], ea_e_dist_lonlat[4], ea_e_dist_lonlat[6]]))
#ea_e_alt1_lons = (ma.concatenate([ea_e_alt1_lonlat[0], ea_e_alt1_lonlat[2], ea_e_alt1_lonlat[4], ea_e_alt1_lonlat[6]]))
#ea_e_alt2_lons = (ma.concatenate([ea_e_alt2_lonlat[0], ea_e_alt2_lonlat[2], ea_e_alt2_lonlat[4], ea_e_alt2_lonlat[6]]))
#ea_e_dist_lats = (ma.concatenate([ea_e_dist_lonlat[1], ea_e_dist_lonlat[3], ea_e_dist_lonlat[5], ea_e_dist_lonlat[7]]))
#ea_e_alt1_lats = (ma.concatenate([ea_e_alt1_lonlat[1], ea_e_alt1_lonlat[3], ea_e_alt1_lonlat[5], ea_e_alt1_lonlat[7]]))
#ea_e_alt2_lats = (ma.concatenate([ea_e_alt2_lonlat[1], ea_e_alt2_lonlat[3], ea_e_alt2_lonlat[5], ea_e_alt2_lonlat[7]]))
#ea_e_lons = ma.concatenate([ea_e_dist_lons, ea_e_alt1_lons, ea_e_alt2_lons])
#ea_e_lats = ma.concatenate([ea_e_dist_lats, ea_e_alt1_lats, ea_e_alt2_lats])#

#na_w_dist_lons = (ma.concatenate([na_w_dist_lonlat[0], na_w_dist_lonlat[2], na_w_dist_lonlat[4], na_w_dist_lonlat[6]]))
#na_w_alt1_lons = (ma.concatenate([na_w_alt1_lonlat[0], na_w_alt1_lonlat[2], na_w_alt1_lonlat[4], na_w_alt1_lonlat[6]]))
#na_w_alt2_lons = (ma.concatenate([na_w_alt2_lonlat[0], na_w_alt2_lonlat[2], na_w_alt2_lonlat[4], na_w_alt2_lonlat[6]]))
#na_w_dist_lats = (ma.concatenate([na_w_dist_lonlat[1], na_w_dist_lonlat[3], na_w_dist_lonlat[5], na_w_dist_lonlat[7]]))
#na_w_alt1_lats = (ma.concatenate([na_w_alt1_lonlat[1], na_w_alt1_lonlat[3], na_w_alt1_lonlat[5], na_w_alt1_lonlat[7]]))
#na_w_alt2_lats = (ma.concatenate([na_w_alt2_lonlat[1], na_w_alt2_lonlat[3], na_w_alt2_lonlat[5], na_w_alt2_lonlat[7]]))
#na_w_lons = ma.concatenate([na_w_dist_lons, na_w_alt1_lons, na_w_alt2_lons])
#na_w_lats = ma.concatenate([na_w_dist_lats, na_w_alt1_lats, na_w_alt2_lats])#

#na_e_dist_lons = (ma.concatenate([na_e_dist_lonlat[0], na_e_dist_lonlat[2], na_e_dist_lonlat[4], na_e_dist_lonlat[6]]))
#na_e_alt1_lons = (ma.concatenate([na_e_alt1_lonlat[0], na_e_alt1_lonlat[2], na_e_alt1_lonlat[4], na_e_alt1_lonlat[6]]))
#na_e_alt2_lons = (ma.concatenate([na_e_alt2_lonlat[0], na_e_alt2_lonlat[2], na_e_alt2_lonlat[4], na_e_alt2_lonlat[6]]))
#na_e_dist_lats = (ma.concatenate([na_e_dist_lonlat[1], na_e_dist_lonlat[3], na_e_dist_lonlat[5], na_e_dist_lonlat[7]]))
#na_e_alt1_lats = (ma.concatenate([na_e_alt1_lonlat[1], na_e_alt1_lonlat[3], na_e_alt1_lonlat[5], na_e_alt1_lonlat[7]]))
#na_e_alt2_lats = (ma.concatenate([na_e_alt2_lonlat[1], na_e_alt2_lonlat[3], na_e_alt2_lonlat[5], na_e_alt2_lonlat[7]]))
#na_e_lons = ma.concatenate([na_e_dist_lons, na_e_alt1_lons, na_e_alt2_lons])
#na_e_lats = ma.concatenate([na_e_dist_lats, na_e_alt1_lats, na_e_alt2_lats])#
#

#commonlons_ea_w, commonlats_ea_w, commonlai_ea_w, commonndvi_ea_w , commontcf_ea_w, commonstate_ea_w = [],[],[],[],[],[]
#lencount_ea_w=0
#for elem in ma.nonzero(ea_w_lons)[0]:
#    commonlons_ea_w.append(data_lons_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])][0])
#    commonlats_ea_w.append(data_lats_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])][0])
#    commonlai_ea_w.append(data_lai_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])][0])
#    commonndvi_ea_w.append(data_ndvi_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])][0])
#    commontcf_ea_w.append(data_tcf_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])][0])
#    commonstate_ea_w.append(data_state_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])][0])#
#
#

#    if len(data_lons_ea_west[(data_lons_ea_west==ea_w_lons[elem]) & (data_lats_ea_west==ea_w_lats[elem])])>1: lencount_ea_w +=1
#print lencount_ea_w#

#laiarr_ea_w = []
#for elem in np.array(commonlai_ea_w)[:]*0.0001: 
#    if (elem>-1.): laiarr_ea_w.append(float("{:10.2f}".format(-1*elem)))
#    else:  laiarr_ea_w.append(float("{:10.2f}".format(elem)))
#ndviarr_ea_w = []
#for elem in np.array(commonndvi_ea_w)[:]*0.0001: 
#    if (elem>-1.): ndviarr_ea_w.append(float("{:10.2f}".format(-1*elem)))
#    else:  ndviarr_ea_w.append(float("{:10.2f}".format(elem)))
#        
#commonlons_ea_e, commonlats_ea_e, commonlai_ea_e, commonndvi_ea_e, commontcf_ea_e, commonstate_ea_e = [],[],[],[],[],[]
#lencount_ea_e=0
#for elem in ma.nonzero(ea_e_lons)[0]:
#    commonlons_ea_e.append(data_lons_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])][0])
#    commonlats_ea_e.append(data_lats_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])][0])
#    commonlai_ea_e.append(data_lai_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])][0])
#    commonndvi_ea_e.append(data_ndvi_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])][0])
#    commontcf_ea_e.append(data_tcf_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])][0])
#    commonstate_ea_e.append(data_state_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])][0])#

#    if len(data_lons_ea_east[(data_lons_ea_east==ea_e_lons[elem]) & (data_lats_ea_east==ea_e_lats[elem])])>1: lencount_ea_e +=1
#print lencount_ea_e#

#laiarr_ea_e = []
#for elem in np.array(commonlai_ea_e)[:]*0.0001: 
#    if (elem>-1.): laiarr_ea_e.append(float("{:10.2f}".format(-1*elem)))
#    else:  laiarr_ea_e.append(float("{:10.2f}".format(elem)))
#ndviarr_ea_e = []
#for elem in np.array(commonndvi_ea_e)[:]*0.0001: 
#    if (elem>-1.): ndviarr_ea_e.append(float("{:10.2f}".format(-1*elem)))
#    else:  ndviarr_ea_e.append(float("{:10.2f}".format(elem)))
#        
#        
#commonlons_na_w, commonlats_na_w, commonlai_na_w, commonndvi_na_w, commontcf_na_w, commonstate_na_w = [],[],[],[],[],[]
#lencount_na_w=0
#for elem in ma.nonzero(na_w_lons)[0]:
#    commonlons_na_w.append(data_lons_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])][0])
#    commonlats_na_w.append(data_lats_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])][0])
#    commonlai_na_w.append(data_lai_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])][0])
#    commonndvi_na_w.append(data_ndvi_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])][0])
#    commontcf_na_w.append(data_tcf_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])][0])
#    commonstate_na_w.append(data_state_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])][0])#

#    if len(data_lons_na_west[(data_lons_na_west==na_w_lons[elem]) & (data_lats_na_west==na_w_lats[elem])])>1: lencount_na_w +=1
#print lencount_na_w#

#laiarr_na_w = []
#for elem in np.array(commonlai_na_w)[:]*0.0001: 
#    if (elem>-1.): laiarr_na_w.append(float("{:10.2f}".format(-1*elem)))
#    else:  laiarr_na_w.append(float("{:10.2f}".format(elem)))
#ndviarr_na_w = []
#for elem in np.array(commonndvi_na_w)[:]*0.0001: 
#    if (elem>-1.): ndviarr_na_w.append(float("{:10.2f}".format(-1*elem)))
#    else:  ndviarr_na_w.append(float("{:10.2f}".format(elem)))
#        
#commonlons_na_e, commonlats_na_e, commonlai_na_e, commonndvi_na_e, commontcf_na_e, commonstate_na_e = [],[],[],[],[],[]
#lencount_na_e=0
#for elem in ma.nonzero(na_e_lons)[0]:
#    commonlons_na_e.append(data_lons_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])][0])
#    commonlats_na_e.append(data_lats_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])][0])
#    commonlai_na_e.append(data_lai_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])][0])
#    commonndvi_na_e.append(data_ndvi_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])][0])
#    commontcf_na_e.append(data_tcf_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])][0])
#    commonstate_na_e.append(data_state_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])][0])#

#    if len(data_lons_na_east[(data_lons_na_east==na_e_lons[elem]) & (data_lats_na_east==na_e_lats[elem])])>1: lencount_na_e +=1
#print lencount_na_e#

#laiarr_na_e = []
#for elem in np.array(commonlai_na_e)[:]*0.0001: 
#    if (elem>-1.): laiarr_na_e.append(float("{:10.2f}".format(-1*elem)))
#    else:  laiarr_na_e.append(float("{:10.2f}".format(elem)))
#ndviarr_na_e = []
#for elem in np.array(commonndvi_na_e)[:]*0.0001: 
#    if (elem>-1.): ndviarr_na_e.append(float("{:10.2f}".format(-1*elem)))
#    else:  ndviarr_na_e.append(float("{:10.2f}".format(elem)))
#        #

#commonlons_bor = np.array(commonlons_ea_e+commonlons_ea_w+commonlons_na_e+commonlons_na_w)
#commonlats_bor = np.array(commonlats_ea_e+commonlats_ea_w+commonlats_na_e+commonlats_na_w)
#commonlai_bor = np.concatenate([laiarr_ea_e, laiarr_ea_w, laiarr_na_e, laiarr_na_w])
#commonndvi_bor = np.concatenate([ndviarr_ea_e, ndviarr_ea_w, ndviarr_na_e, ndviarr_na_w])
#commontcf_bor = np.array(commontcf_ea_e+commontcf_ea_w+commontcf_na_e+commontcf_na_w)
#commonstate_bor = np.array(commonstate_ea_e+commonstate_ea_w+commonstate_na_e+commonstate_na_w)
#print len(commonlai_bor), len(commonlons_bor), len(commonlats_bor), len(commonndvi_bor), len(commontcf_bor), len(commonstate_bor)
##np.savetxt('commonlons_bor.txt', commonlons_bor)
##np.savetxt('commonlats_bor.txt', commonlats_bor)
##np.savetxt('commonlai_bor.txt', commonlai_bor)
##np.savetxt('commonndvi_bor.txt', commonndvi_bor)
##np.savetxt('commontcf_bor.txt', commontcf_bor)
##np.savetxt('commonstate_bor.txt', commonstate_bor)

