{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap, shiftgrid\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import datetime\n",
    "import rpy2.robjects as robjects\n",
    "import rpy2.robjects.packages as rpackages # import rpy2's package module\n",
    "from rpy2.robjects.packages import importr # import R's utility package\n",
    "import netCDF4\n",
    "import sys\n",
    "import math \n",
    "import numpy.ma as ma\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# sys.path.append(\"/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/\") # Edited .zshrc profile with PYTHONPATH \n",
    "from Functions import plot_north_polar\n",
    "from pandas import DataFrame \n",
    "from collections import Counter\n",
    "\n",
    "\n",
    "r = robjects.r\n",
    "# Some R packages\n",
    "mgcv = importr('mgcv')\n",
    "importr(\"data.table\")\n",
    "utils = importr('utils')\n",
    "# select a mirror for R packages\n",
    "utils.chooseCRANmirror(ind=1) # select the first mirror in the list\n",
    "from rpy2.robjects import FloatVector\n",
    "import rpy2.robjects.numpy2ri # For the conversion of numpy objects into R \n",
    "rpy2.robjects.numpy2ri.activate()\n",
    "silvermantest = importr('silvermantest')\n",
    "\n",
    "\n",
    "sns.set_style(style=\"white\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#final arrays with lon, lat of all alternative states \n",
    "commonlons_bor = np.loadtxt('commonlons_bor.txt')\n",
    "commonlats_bor = np.loadtxt('commonlats_bor.txt')\n",
    "commonlai_bor = np.loadtxt('commonlai_bor.txt')\n",
    "commonndvi_bor = np.loadtxt('commonndvi_bor.txt')\n",
    "commontcf_bor = np.loadtxt('commontcf_bor.txt')\n",
    "commonstate_bor = np.loadtxt('commonstate_bor.txt')\n",
    "orderednums = np.loadtxt('orderednums.txt')\n",
    "\n",
    "\n",
    "# let's divide TCF into the three states\n",
    "tcf_less = [elem for elem in commontcf_bor if elem<20]\n",
    "tcf_open = [elem for elem in commontcf_bor if (elem>=20 and elem<45)]\n",
    "tcf_fore = [elem for elem in commontcf_bor if elem>=45]\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "colbare = '#332288'\n",
    "colsavab = '#88CCEE'\n",
    "colsava = '#44AA99'\n",
    "colfore = '#117733'\n",
    "\n",
    "alt_colbare = '#999933'\n",
    "alt_colsavab = '#CC6677' \n",
    "alt_colsavaf = '#AA4499'\n",
    "alt_colfore = '#882255'\n",
    "\n",
    "\n",
    "def setcolor(x, color): # for basemap plots-> color parallels and meridians\n",
    "    for m in x:\n",
    "         for t in x[m][1]:\n",
    "             t.set_color(color)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# execute this only if you want to re compute the variables\n",
    "\n",
    "path = '/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Boreal_New_Classification/'\n",
    "\n",
    "na_e_dist_lonlat = np.load(path+'na_e_dist_lolat.npy')\n",
    "na_e_alt1_lonlat = np.load(path+'na_e_alt1_lolat.npy')\n",
    "na_e_alt2_lonlat = np.load(path+'na_e_alt2_lolat.npy')\n",
    "na_w_dist_lonlat = np.load(path+'na_w_dist_lolat.npy')\n",
    "na_w_alt1_lonlat = np.load(path+'na_w_alt1_lolat.npy')\n",
    "na_w_alt2_lonlat = np.load(path+'na_w_alt2_lolat.npy')\n",
    "ea_w_dist_lonlat = np.load(path+'ea_w_dist_lolat.npy')\n",
    "ea_w_alt1_lonlat = np.load(path+'ea_w_alt1_lolat.npy')\n",
    "ea_w_alt2_lonlat = np.load(path+'ea_w_alt2_lolat.npy')\n",
    "ea_e_dist_lonlat = np.load(path+'ea_e_dist_lolat.npy')\n",
    "ea_e_alt1_lonlat = np.load(path+'ea_e_alt1_lolat.npy')\n",
    "ea_e_alt2_lonlat = np.load(path+'ea_e_alt2_lolat.npy')\n",
    "\n",
    "\n",
    "# concatenate lonlat arrays of fire disturbed and equivalent tree cover arrays\n",
    "# create same shape arrays filled with increasing numbers\n",
    "altlist = [na_e_dist_lonlat, na_w_dist_lonlat, ea_e_dist_lonlat, ea_w_dist_lonlat]\n",
    "altlon = np.concatenate([np.concatenate([ma.compressed(elem[i*2]) for i in range(4)]) for elem in altlist])\n",
    "altlat = np.concatenate([np.concatenate([ma.compressed(elem[i*2+1]) for i in range(4)]) for elem in altlist])\n",
    "altnum = np.concatenate([np.concatenate([np.ones(ma.compressed(elem[i*2+1]).shape)*i for i in range(4)]) for elem in altlist])\n",
    "\n",
    "fdlist = [na_e_alt1_lonlat, na_w_alt1_lonlat, ea_e_alt1_lonlat, ea_w_alt1_lonlat, \\\n",
    "           na_e_alt2_lonlat, na_w_alt2_lonlat, ea_e_alt2_lonlat, ea_w_alt2_lonlat] \n",
    "fdlon = np.concatenate([np.concatenate([ma.compressed(elem[i*2]) for i in range(4)]) for elem in fdlist])\n",
    "fdlat = np.concatenate([np.concatenate([ma.compressed(elem[i*2+1]) for i in range(4)]) for elem in fdlist])\n",
    "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])\n",
    "\n",
    "# concatenate into single arrays to have everything together\n",
    "altfdlon = np.concatenate([altlon,fdlon])\n",
    "altfdlat = np.concatenate([altlat,fdlat])\n",
    "altfdnum = np.concatenate([altnum,fdnum])\n",
    "\n",
    "#create cmap for alternative states\n",
    "colseq = [colbare,colsavab,colsava,colfore,alt_colbare,alt_colsavab,alt_colsavaf,alt_colfore]\n",
    "current_palette = sns.color_palette(colseq)\n",
    "altfdmap = ListedColormap(current_palette)\n",
    "\n",
    "# sort a lons of the altfd array so that they match the commonlons_bor, and same for the numerical array necessary for colouring\n",
    "objmap = {altfdlon[i]: altfdnum[i] for i in range(len(altfdlon))}\n",
    "orderednums = [objmap[id] for id in commonlons_bor]\n",
    "#np.savetxt('orderednums.txt', orderednums)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let's check bimodality of every combination (open treeless, open forest), for every region->next put the regions together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "minimum peaks distance 18.1917808219\n",
      "all peaks distances [21.956947162426612, 18.19178082191781]\n",
      "Number of gridcells with variability >18.1917808219 1\n",
      "Mean standard deviation 5.77030467449 \n",
      "Median standard deviation 5.33760245803 \n",
      "Standard deviation of standard deviation 2.6953706488\n"
     ]
    }
   ],
   "source": [
    "# creates the variables for comparisons\n",
    "import pandas as pd\n",
    "zone = '_ea_east'\n",
    "path_to_pickles = '/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/'\n",
    "bare_classified_ea_e = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')\n",
    "savanna_classified_ea_e = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')\n",
    "forest_classified_ea_e = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_less_common_ea_e = np.concatenate([treeless_open_common_ea_e, open_treeless_common_ea_e])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_fore_common_ea_e = np.concatenate([forest_open_common_ea_e, open_forest_common_ea_e])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "#r.assign('tcf_bor_op_less_ea_e', np.array(bare_open_tcf_comparison_ea_e))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_ea_e),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "#r.assign('tcf_bor_op_fore_ea_e', np.array(forest_open_tcf_comparison_ea_e))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_ea_e),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "#EA_W\n",
    "zone = '_ea_west_russia'\n",
    "bare_classified_ea_w = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')\n",
    "savanna_classified_ea_w = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')\n",
    "forest_classified_ea_w = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_less_common_ea_w = np.concatenate([treeless_open_common_ea_w, open_treeless_common_ea_w])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_fore_common_ea_w = np.concatenate([forest_open_common_ea_w, open_forest_common_ea_w])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "#r.assign('tcf_bor_op_less_ea_w', np.array(bare_open_tcf_comparison_ea_w))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_ea_w),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "#r.assign('tcf_bor_op_fore_ea_w', np.array(forest_open_tcf_comparison_ea_w))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_ea_w),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "\n",
    "# NA_E\n",
    "zone = '_na_east'\n",
    "bare_classified_na_e = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')\n",
    "savanna_classified_na_e = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')\n",
    "forest_classified_na_e = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_less_common_na_e = np.concatenate([treeless_open_common_na_e, open_treeless_common_na_e])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_fore_common_na_e = np.concatenate([forest_open_common_na_e, open_forest_common_na_e])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "#r.assign('tcf_bor_op_less_na_e', np.array(bare_open_tcf_comparison_na_e))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_na_e),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "#r.assign('tcf_bor_op_fore_na_e', np.array(forest_open_tcf_comparison_na_e))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_na_e),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "#NA_W\n",
    "zone = '_na_west'\n",
    "bare_classified_na_w = pd.read_pickle(path_to_pickles+'bare'+zone+'_classified_with_SoilT.pkl')\n",
    "savanna_classified_na_w = pd.read_pickle(path_to_pickles+'savanna'+zone+'_classified_with_SoilT.pkl')\n",
    "forest_classified_na_w = pd.read_pickle(path_to_pickles+'forest'+zone+'_classified_with_SoilT.pkl')\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_less_common_na_w = np.concatenate([treeless_open_common_na_w, open_treeless_common_na_w])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "\n",
    "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])])\n",
    "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])])\n",
    "\n",
    "op_fore_common_na_w = np.concatenate([forest_open_common_na_w, open_forest_common_na_w])\n",
    "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]\n",
    "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)])\n",
    "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)])\n",
    "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)])\n",
    "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))])\n",
    "\n",
    "# not enough points to compute\n",
    "#r.assign('tcf_bor_op_less_na_w', np.array(bare_open_tcf_comparison_na_w))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_less_na_w),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "#r.assign('tcf_bor_op_fore_na_w', np.array(forest_open_tcf_comparison_na_w))\n",
    "#r('stest <- silverman.test(as.numeric(tcf_bor_op_fore_na_w),1, adjust=TRUE)');\n",
    "#print r('stest')\n",
    "\n",
    "\n",
    "\n",
    "# lon lat for boreal open treeless and open forest\n",
    "\n",
    "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])\n",
    "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])\n",
    "\n",
    "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])\n",
    "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])\n",
    "\n",
    "\n",
    "#entire boreal area\n",
    "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])\n",
    "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])\n",
    "\n",
    "#open_treeless_bor = np.loadtxt(\"open_treless_bor.txt\")\n",
    "#forest_open_bor = np.loadtxt(\"forest_open_bor.txt\")\n",
    "# let's divide them into the single states\n",
    "\n",
    "alt_less = [elem for elem in open_treeless_bor if elem<20]\n",
    "alt_open_l = [elem for elem in open_treeless_bor if (elem>=20)]\n",
    "alt_open_f = [elem for elem in forest_open_bor if (elem<45)]\n",
    "alt_fore = [elem for elem in forest_open_bor if elem>=45]\n",
    "\n",
    "# Fitting of the histograms with KDE\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.distributions.mixture_rvs import mixture_rvs\n",
    "from statsmodels.nonparametric.kernel_density import KDEMultivariate as kdem\n",
    "import peakutils\n",
    "kde = sm.nonparametric.KDEUnivariate(open_treeless_bor)\n",
    "kde.fit(bw=4)\n",
    "kde2 = sm.nonparametric.KDEUnivariate(forest_open_bor)\n",
    "kde2.fit(bw=4)\n",
    "\n",
    "indexes = peakutils.indexes(kde.density, thres=0.02)#, thres=0.02/max(cb), min_dist=100)\n",
    "indexes2 = peakutils.indexes(kde2.density, thres=0.02)#, thres=0.02/max(cb), min_dist=100)\n",
    "peaks_x = [kde.support[i] for i in indexes] \n",
    "peaks_y = [kde.density[i] for i in indexes]\n",
    "peaks_x2 = [kde2.support[i] for i in indexes2] \n",
    "peaks_y2 = [kde2.density[i] for i in indexes2]\n",
    "\n",
    "\n",
    "\n",
    "min_peak_dist = np.min([peaks_x[1]-peaks_x[0], peaks_x2[1]-peaks_x2[0]])\n",
    "print \"minimum peaks distance\", min_peak_dist\n",
    "print \"all peaks distances\", [peaks_x[1]-peaks_x[0], peaks_x2[1]-peaks_x2[0]]\n",
    "\n",
    "\n",
    "\n",
    "#Variability\n",
    "\n",
    "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') \n",
    "\n",
    "tcf_years = ma.getdata(nc_tcf_years.variables['vcff'][:])\n",
    "tcf_y_lons, tcf_y_lats = nc_tcf_years.variables['lon'][:], nc_tcf_years.variables['lat'][:]\n",
    "\n",
    "tcf_y_bor = []\n",
    "for elem in tcf_years:\n",
    "    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))]))\n",
    "tcf_y_arr = np.array(tcf_y_bor)\n",
    "\n",
    "\n",
    "max_array, min_array = [], []\n",
    "for i in range(tcf_y_arr.shape[1]):\n",
    "    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])\n",
    "    min_array.append(np.min(tcf_y_arr[:,i]))\n",
    "    \n",
    "    \n",
    "vari_array = (np.array(max_array)-np.array(min_array))\n",
    "np.mean(vari_array)\n",
    "\n",
    "std_var_arr = np.array([np.std(tcf_y_arr[:,i]) for i in range(tcf_y_arr.shape[1])])\n",
    "print \"Number of gridcells with variability >\"+str(min_peak_dist), len([elem for elem in std_var_arr if elem>min_peak_dist])\n",
    "import statistics as st\n",
    "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\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Silvermantest: Testing the hypothesis if the number of modes is <=  1 \n",
      "The P-Value is  0.003845 \n",
      "\n",
      "Silvermantest: Testing the hypothesis if the number of modes is <=  1 \n",
      "The P-Value is  0.021033 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Bimodality tests using silverman's test for the Open woodland + Treeless states \n",
    "r.assign('open_treeless_bor', np.array(open_treeless_bor))\n",
    "r('stest <- silverman.test(as.numeric(open_treeless_bor),1, adjust=TRUE)');\n",
    "print r('stest')\n",
    "\n",
    "# Bimodality tests using silverman's test for the Open woodland + Forest states \n",
    "r.assign('forest_open_bor', np.array(forest_open_bor))\n",
    "r('stest <- silverman.test(as.numeric(forest_open_bor),1, adjust=TRUE)');\n",
    "print r('stest')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Silvermantest: Testing the hypothesis if the number of modes is <=  2 \n",
      "The P-Value is  0.003003003 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Silverman's test for trimodality of alternative states \n",
    "r.assign('tcf_bor', np.array(commontcf_bor))\n",
    "r('stest <- silverman.test(as.numeric(tcf_bor),2, adjust=FALSE)');\n",
    "print r('stest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Some plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Bimodality of alternative states categories: Open woodland + Treeless and Open woodland + Forest\n",
    "# Histograms and silverman's test for bimodality\n",
    "fig = plt.figure(figsize=(20,9))\n",
    "ax = fig.add_subplot(121)\n",
    "ax.set_title('Alternative Treeless / Open Woodland States', fontsize=16)\n",
    "hi = plt.hist([alt_less, alt_open_l], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#7570b3', '#d95f02']\\\n",
    "        , label=[\"Treeless\", \"Open woodland\"])\n",
    "sns.despine()\n",
    "plt.xlim(0,50)\n",
    "plt.xlabel(\"Tree Cover Fraction - [%]\", fontsize=14)\n",
    "plt.ylabel(\"Normalised Frequency\", fontsize=14)\n",
    "plt.legend(fontsize=14,ncol=3, loc='upper center')\n",
    "ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\\leq$'+' 1 \\n The P-Value is '+r'$<$'+' 0.001', \\\n",
    "            xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "ax.set_title('Alternative Forest / Open Woodland States', fontsize=16)\n",
    "hi = plt.hist([alt_open_f, alt_fore], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#d95f02', '#1b9e77']\\\n",
    "        , label=[\"Open woodland\", \"Forest\"])\n",
    "sns.despine()\n",
    "plt.xlim(20,80)\n",
    "plt.ylim(0,0.12)\n",
    "plt.xlabel(\"Tree Cover Fraction - [%]\", fontsize=14)\n",
    "plt.ylabel(\"Normalised Frequency\", fontsize=14)\n",
    "plt.legend(fontsize=14,ncol=3, loc='upper center')\n",
    "ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\\leq$'+' 1 \\n The P-Value is '+r'$<$'+' 0.01 ', \\\n",
    "            xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)\n",
    "\n",
    "#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Alternative_States_TCF_Bimodality.pdf', dpi=600, transparent = True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fitting of the distributions of the two classes of alternative states using KDE, \n",
    "# peaks distance and variability comparison\n",
    "\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title(\"Alternative States TCF KDE\", fontsize=16)\n",
    "ax.hist(open_treeless_bor, bins=10, normed=True, color='#7570b3', alpha=0.7, label = 'Treeless / Open woodland States')\n",
    "ax.plot(kde.support, kde.density, lw=2, color='black', label = 'KDE estimates');\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(111)\n",
    "ax.hist(forest_open_bor, bins=10, normed=True, color='#1b9e77', alpha=0.7, label='Forest / Open woodland States')\n",
    "ax.plot(kde2.support, kde2.density, lw=2, color='black');\n",
    "plt.xlabel(\"Tree Cover Fraction - [%]\", fontsize=14)\n",
    "plt.ylabel(\"Normalised Frequency\", fontsize=14)\n",
    "sns.despine()\n",
    "plt.xlim(-10,90)\n",
    "\n",
    "\n",
    "\n",
    "plt.plot(peaks_x, peaks_y, 'o', color='#d95f02', label = 'KDE peaks')\n",
    "plt.plot(peaks_x2, peaks_y2, 'o', color='#d95f02')\n",
    "\n",
    "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",
    "            '\\n Minimum distance between KDE peaks is '+str(\"{:1.2f}\".format(min_peak_dist)), \\\n",
    "            xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)\n",
    "plt.legend(fontsize=14)\n",
    "#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Alternative_States_TCF_Variability.pdf', dpi=600, transparent = True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Bimodality of alternative states, with Silverman's and parallele to TCF plots\n",
    "fig = plt.figure(figsize=(15,11))\n",
    "ax = fig.add_subplot(221)\n",
    "ax.annotate(\"Open Woodland / Treeless TCF\",  xy=(0.5,1.1), xycoords='axes fraction', ha='center', fontsize=18)\n",
    "m1 = Basemap(projection='npaeqd',boundinglat=45,lon_0=-90,resolution='l', round=True)\n",
    "\n",
    "# map extra features # labels are at intersection the left, right, top or bottom of the plot\n",
    "par = m1.drawparallels([45.,55.,65.,75.,90.],labels=[False,False,False,False], color='black', fontsize=16)\n",
    "mer = m1.drawmeridians(np.arange(-135.,179.,45.),labels=[False,False,True,True], color='black', fontsize=16)\n",
    "m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'\n",
    "m1.drawcoastlines(linewidth=1.5)\n",
    "m1.fillcontinents(color='white',lake_color='black')\n",
    "m1_sub_lons, m1_sub_lats = m1(open_treeless_bor_lon,open_treeless_bor_lat)\n",
    "cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=open_treeless_bor, cmap = 'viridis',\\\n",
    "                 alpha=0.8, s=20, zorder=3)\n",
    "m1.colorbar(cbs,\"right\", ticks=np.arange(0,90,10),size=\"5%\", pad=\"15%\")\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(222)\n",
    "hi = plt.hist([alt_less, alt_open_l], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#7570b3', '#d95f02']\\\n",
    "        , label=[\"Treeless\", \"Open woodland\"])\n",
    "sns.despine()\n",
    "plt.xlim(0,50)\n",
    "plt.xlabel(\"Tree Cover Fraction - [%]\", fontsize=14)\n",
    "plt.ylabel(\"Normalised Frequency\", fontsize=14)\n",
    "plt.legend(fontsize=14,ncol=3, loc='upper center')\n",
    "ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\\leq$'+' 1 \\n The P-Value is  0.001', \\\n",
    "            xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(223)\n",
    "ax.annotate(\"Forest / Open Woodland TCF\",  xy=(0.5,1.1), xycoords='axes fraction', ha='center', fontsize=18)\n",
    "m1 = Basemap(projection='npaeqd',boundinglat=45,lon_0=-90,resolution='l', round=True)\n",
    "\n",
    "# map extra features # labels are at intersection the left, right, top or bottom of the plot\n",
    "par = m1.drawparallels([45.,55.,65.,75.,90.],labels=[False,False,False,False], color='black', fontsize=16)\n",
    "mer = m1.drawmeridians(np.arange(-135.,179.,45.),labels=[False,False,True,True], color='black', fontsize=16)\n",
    "m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'\n",
    "m1.drawcoastlines(linewidth=1.5)\n",
    "m1.fillcontinents(color='white',lake_color='black')\n",
    "m1_sub_lons, m1_sub_lats = m1(forest_open_bor_lon,forest_open_bor_lat)\n",
    "cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=forest_open_bor, cmap = 'viridis',\\\n",
    "                 alpha=0.8, s=20, zorder=3)\n",
    "m1.colorbar(cbs,\"right\", ticks=np.arange(0,90,10),size=\"5%\", pad=\"15%\")\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(224)\n",
    "hi = plt.hist([alt_open_f, alt_fore], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#d95f02', '#1b9e77']\\\n",
    "        , label=[\"Open woodland\", \"Forest\"])\n",
    "sns.despine()\n",
    "plt.xlim(20,80)\n",
    "plt.ylim(0,0.12)\n",
    "plt.xlabel(\"Tree Cover Fraction - [%]\", fontsize=14)\n",
    "plt.ylabel(\"Normalised Frequency\", fontsize=14)\n",
    "plt.legend(fontsize=14,ncol=3, loc='upper center')\n",
    "ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\\leq$'+' 1 \\n The P-Value is 0.01 ', \\\n",
    "            xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TCF plot and histogram for trimodality with silverman's test\n",
    "fig = plt.figure(figsize=(20,9))\n",
    "ax = fig.add_subplot(121)\n",
    "ax.annotate(\"Alternative States TCF\",  xy=(0.5,1.1), xycoords='axes fraction', ha='center', fontsize=18)\n",
    "m1 = Basemap(projection='npaeqd',boundinglat=45,lon_0=-90,resolution='l', round=True)\n",
    "\n",
    "# map extra features # labels are at intersection the left, right, top or bottom of the plot\n",
    "par = m1.drawparallels([45.,55.,65.,75.,90.],labels=[False,False,False,False], color='black', fontsize=16)\n",
    "mer = m1.drawmeridians(np.arange(-135.,179.,45.),labels=[False,False,True,True], color='black', fontsize=16)\n",
    "\n",
    "m1.drawmapboundary(fill_color='#02c5e9')#black'#02d8e9'\n",
    "m1.drawcoastlines(linewidth=1.5)\n",
    "m1.fillcontinents(color='white',lake_color='black')\n",
    "\n",
    "# set the color of the parallel and meridian lines\n",
    "#setcolor(par,'black')\n",
    "#setcolor(mer,'black')\n",
    "\n",
    "\n",
    "m1_sub_lons, m1_sub_lats = m1(commonlons_bor,commonlats_bor)\n",
    "cbs = m1.scatter(m1_sub_lons,m1_sub_lats, marker='h', c=commontcf_bor, cmap = 'viridis',\\\n",
    "                 alpha=0.8, s=20, zorder=3)\n",
    "\n",
    "m1.colorbar(cbs,\"right\", ticks=np.arange(0,90,10),size=\"5%\", pad=\"15%\")\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "ax.set_title('Alternative States TCF Multimodality Test', fontsize=18)\n",
    "hi = plt.hist([tcf_less, tcf_open, tcf_fore], bins=20, rwidth=2, alpha=0.7, normed=True, color = ['#7570b3', '#d95f02', '#1b9e77']\\\n",
    "        , label=[\"Treeless\", \"Open woodland\", \"Forest\"])\n",
    "sns.despine()\n",
    "plt.xlim(0,80)\n",
    "plt.xlabel(\"Tree Cover Fraction - [%]\", fontsize=14)\n",
    "plt.ylabel(\"Normalised Frequency\", fontsize=14)\n",
    "plt.legend(fontsize=14,ncol=3, loc='upper center')\n",
    "ax.annotate('Silvermantest: Testing the hypothesis if the number of modes is '+r'$\\leq$'+' 2 \\n The P-Value is  0.001001001 ', \\\n",
    "            xy=(0.5,0.7), xycoords='axes fraction', ha='center', fontsize=14)\n",
    "\n",
    "#plt.savefig('/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Pic/Alternative_States_TCF_Multimodality.pdf', dpi=600, transparent = True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# results of silverman's test with different number of modes\n",
    "stests_ol = [0.002005, 0.3503504, 0.02102102, 0.4614615, 0.4254254] # results pf silvermantests with 1, 2, 3, 4, 5 modes\n",
    "stests_of = [0.012927, 0.7897898, 0.4744745, 0.4764765, 0.6216216]\n",
    "x_stests = [1,2,3,4,5]\n",
    "null_line = np.arange(1,5, 0.01)\n",
    "x_line = 0.05*np.ones(null_line.shape)\n",
    "fig = plt.figure()\n",
    "plt.plot(null_line, x_line, 'r-')\n",
    "plt.plot(x_stests, stests_ol, '-o', color = '#7570b3', label = 'Treeless / Open woodland')\n",
    "plt.plot(x_stests, stests_of, '-o', color = '#1b9e77', label = 'Forest / Open woodland')\n",
    "\n",
    "plt.xlim(0.8,5.2)\n",
    "plt.ylim(-0.02,1.0)\n",
    "plt.xlabel('Number of modes in null hypothesis of silvermantest')\n",
    "plt.ylabel('p-value')\n",
    "plt.legend()\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# data load execute to re-create the initial common_variables\n",
    "\n",
    "#dataforr_ea_west = np.loadtxt(\"boreal_ea_west_database_lai_ndvi_included.txt\")\n",
    "#data_lons_ea_west = dataforr_ea_west[:,-2]\n",
    "#data_lats_ea_west = dataforr_ea_west[:,-1]\n",
    "#data_lai_ea_west = dataforr_ea_west[:,10]\n",
    "#data_ndvi_ea_west = dataforr_ea_west[:,11]\n",
    "#data_tcf_ea_west = dataforr_ea_west[:,0]\n",
    "#data_state_ea_west = dataforr_ea_west[:,-6]\n",
    "\n",
    "\n",
    "\n",
    "#dataforr_ea_east = np.loadtxt(\"boreal_ea_east_database_lai_ndvi_included.txt\")\n",
    "#data_lons_ea_east = dataforr_ea_east[:,-2]\n",
    "#data_lats_ea_east = dataforr_ea_east[:,-1]\n",
    "#data_lai_ea_east = dataforr_ea_east[:,10]\n",
    "#data_ndvi_ea_east = dataforr_ea_east[:,11]\n",
    "#data_tcf_ea_east = dataforr_ea_east[:,0]\n",
    "#data_state_ea_east = dataforr_ea_east[:,-6]\n",
    "\n",
    "\n",
    "#dataforr_na_east = np.loadtxt(\"boreal_na_east_database_lai_ndvi_included.txt\")\n",
    "#data_lons_na_east = dataforr_na_east[:,-2]\n",
    "#data_lats_na_east = dataforr_na_east[:,-1]\n",
    "#data_lai_na_east = dataforr_na_east[:,10]\n",
    "#data_ndvi_na_east = dataforr_na_east[:,11]\n",
    "#data_tcf_na_east = dataforr_na_east[:,0]\n",
    "#data_state_na_east = dataforr_na_east[:,-6]#\n",
    "\n",
    "#dataforr_na_west = np.loadtxt(\"boreal_na_west_database_lai_ndvi_included.txt\")\n",
    "#data_lons_na_west = dataforr_na_west[:,-2]\n",
    "#data_lats_na_west = dataforr_na_west[:,-1]\n",
    "#data_lai_na_west = dataforr_na_west[:,10]\n",
    "#data_ndvi_na_west = dataforr_na_west[:,11]\n",
    "#data_tcf_na_west = dataforr_na_west[:,0]\n",
    "#data_state_na_west = dataforr_na_west[:,-6]#\n",
    "#\n",
    "#\n",
    "\n",
    "## alternative states data#\n",
    "#\n",
    "\n",
    "#path = '/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Boreal_New_Classification/'#\n",
    "\n",
    "#na_e_dist_lonlat = np.load(path+'na_e_dist_lolat.npy')\n",
    "#na_e_alt1_lonlat = np.load(path+'na_e_alt1_lolat.npy')\n",
    "#na_e_alt2_lonlat = np.load(path+'na_e_alt2_lolat.npy')\n",
    "#na_w_dist_lonlat = np.load(path+'na_w_dist_lolat.npy')\n",
    "#na_w_alt1_lonlat = np.load(path+'na_w_alt1_lolat.npy')\n",
    "#na_w_alt2_lonlat = np.load(path+'na_w_alt2_lolat.npy')\n",
    "#ea_w_dist_lonlat = np.load(path+'ea_w_dist_lolat.npy')\n",
    "#ea_w_alt1_lonlat = np.load(path+'ea_w_alt1_lolat.npy')\n",
    "#ea_w_alt2_lonlat = np.load(path+'ea_w_alt2_lolat.npy')\n",
    "#ea_e_dist_lonlat = np.load(path+'ea_e_dist_lolat.npy')\n",
    "#ea_e_alt1_lonlat = np.load(path+'ea_e_alt1_lolat.npy')\n",
    "#ea_e_alt2_lonlat = np.load(path+'ea_e_alt2_lolat.npy')#\n",
    "\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#ea_w_lons = ma.concatenate([ea_w_dist_lons, ea_w_alt1_lons, ea_w_alt2_lons])\n",
    "#ea_w_lats = ma.concatenate([ea_w_dist_lats, ea_w_alt1_lats, ea_w_alt2_lats])#\n",
    "\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#ea_e_lons = ma.concatenate([ea_e_dist_lons, ea_e_alt1_lons, ea_e_alt2_lons])\n",
    "#ea_e_lats = ma.concatenate([ea_e_dist_lats, ea_e_alt1_lats, ea_e_alt2_lats])#\n",
    "\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#na_w_lons = ma.concatenate([na_w_dist_lons, na_w_alt1_lons, na_w_alt2_lons])\n",
    "#na_w_lats = ma.concatenate([na_w_dist_lats, na_w_alt1_lats, na_w_alt2_lats])#\n",
    "\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#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]]))\n",
    "#na_e_lons = ma.concatenate([na_e_dist_lons, na_e_alt1_lons, na_e_alt2_lons])\n",
    "#na_e_lats = ma.concatenate([na_e_dist_lats, na_e_alt1_lats, na_e_alt2_lats])#\n",
    "#\n",
    "\n",
    "#commonlons_ea_w, commonlats_ea_w, commonlai_ea_w, commonndvi_ea_w , commontcf_ea_w, commonstate_ea_w = [],[],[],[],[],[]\n",
    "#lencount_ea_w=0\n",
    "#for elem in ma.nonzero(ea_w_lons)[0]:\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])#\n",
    "#\n",
    "#\n",
    "\n",
    "#    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\n",
    "#print lencount_ea_w#\n",
    "\n",
    "#laiarr_ea_w = []\n",
    "#for elem in np.array(commonlai_ea_w)[:]*0.0001: \n",
    "#    if (elem>-1.): laiarr_ea_w.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  laiarr_ea_w.append(float(\"{:10.2f}\".format(elem)))\n",
    "#ndviarr_ea_w = []\n",
    "#for elem in np.array(commonndvi_ea_w)[:]*0.0001: \n",
    "#    if (elem>-1.): ndviarr_ea_w.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  ndviarr_ea_w.append(float(\"{:10.2f}\".format(elem)))\n",
    "#        \n",
    "#commonlons_ea_e, commonlats_ea_e, commonlai_ea_e, commonndvi_ea_e, commontcf_ea_e, commonstate_ea_e = [],[],[],[],[],[]\n",
    "#lencount_ea_e=0\n",
    "#for elem in ma.nonzero(ea_e_lons)[0]:\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])#\n",
    "\n",
    "#    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\n",
    "#print lencount_ea_e#\n",
    "\n",
    "#laiarr_ea_e = []\n",
    "#for elem in np.array(commonlai_ea_e)[:]*0.0001: \n",
    "#    if (elem>-1.): laiarr_ea_e.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  laiarr_ea_e.append(float(\"{:10.2f}\".format(elem)))\n",
    "#ndviarr_ea_e = []\n",
    "#for elem in np.array(commonndvi_ea_e)[:]*0.0001: \n",
    "#    if (elem>-1.): ndviarr_ea_e.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  ndviarr_ea_e.append(float(\"{:10.2f}\".format(elem)))\n",
    "#        \n",
    "#        \n",
    "#commonlons_na_w, commonlats_na_w, commonlai_na_w, commonndvi_na_w, commontcf_na_w, commonstate_na_w = [],[],[],[],[],[]\n",
    "#lencount_na_w=0\n",
    "#for elem in ma.nonzero(na_w_lons)[0]:\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])#\n",
    "\n",
    "#    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\n",
    "#print lencount_na_w#\n",
    "\n",
    "#laiarr_na_w = []\n",
    "#for elem in np.array(commonlai_na_w)[:]*0.0001: \n",
    "#    if (elem>-1.): laiarr_na_w.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  laiarr_na_w.append(float(\"{:10.2f}\".format(elem)))\n",
    "#ndviarr_na_w = []\n",
    "#for elem in np.array(commonndvi_na_w)[:]*0.0001: \n",
    "#    if (elem>-1.): ndviarr_na_w.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  ndviarr_na_w.append(float(\"{:10.2f}\".format(elem)))\n",
    "#        \n",
    "#commonlons_na_e, commonlats_na_e, commonlai_na_e, commonndvi_na_e, commontcf_na_e, commonstate_na_e = [],[],[],[],[],[]\n",
    "#lencount_na_e=0\n",
    "#for elem in ma.nonzero(na_e_lons)[0]:\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])\n",
    "#    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])#\n",
    "\n",
    "#    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\n",
    "#print lencount_na_e#\n",
    "\n",
    "#laiarr_na_e = []\n",
    "#for elem in np.array(commonlai_na_e)[:]*0.0001: \n",
    "#    if (elem>-1.): laiarr_na_e.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  laiarr_na_e.append(float(\"{:10.2f}\".format(elem)))\n",
    "#ndviarr_na_e = []\n",
    "#for elem in np.array(commonndvi_na_e)[:]*0.0001: \n",
    "#    if (elem>-1.): ndviarr_na_e.append(float(\"{:10.2f}\".format(-1*elem)))\n",
    "#    else:  ndviarr_na_e.append(float(\"{:10.2f}\".format(elem)))\n",
    "#        #\n",
    "\n",
    "#commonlons_bor = np.array(commonlons_ea_e+commonlons_ea_w+commonlons_na_e+commonlons_na_w)\n",
    "#commonlats_bor = np.array(commonlats_ea_e+commonlats_ea_w+commonlats_na_e+commonlats_na_w)\n",
    "#commonlai_bor = np.concatenate([laiarr_ea_e, laiarr_ea_w, laiarr_na_e, laiarr_na_w])\n",
    "#commonndvi_bor = np.concatenate([ndviarr_ea_e, ndviarr_ea_w, ndviarr_na_e, ndviarr_na_w])\n",
    "#commontcf_bor = np.array(commontcf_ea_e+commontcf_ea_w+commontcf_na_e+commontcf_na_w)\n",
    "#commonstate_bor = np.array(commonstate_ea_e+commonstate_ea_w+commonstate_na_e+commonstate_na_w)\n",
    "#print len(commonlai_bor), len(commonlons_bor), len(commonlats_bor), len(commonndvi_bor), len(commontcf_bor), len(commonstate_bor)\n",
    "##np.savetxt('commonlons_bor.txt', commonlons_bor)\n",
    "##np.savetxt('commonlats_bor.txt', commonlats_bor)\n",
    "##np.savetxt('commonlai_bor.txt', commonlai_bor)\n",
    "##np.savetxt('commonndvi_bor.txt', commonndvi_bor)\n",
    "##np.savetxt('commontcf_bor.txt', commontcf_bor)\n",
    "##np.savetxt('commonstate_bor.txt', commonstate_bor)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
