{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import datetime\n",
    "import netCDF4\n",
    "import sys\n",
    "import math \n",
    "\n",
    "# sys.path.append(\"/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/\") # Edited .zshrc profile with PYTHONPATH \n",
    "from pandas import DataFrame \n",
    "from collections import Counter\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# path_to_dataset = /Users/Beniamino/Documents/Datasets/Boreal\n",
    "# execute this only if you want to do the basemap plot\n",
    "\n",
    "#nc005ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.05deg__UHAM-ICDC.nc')\n",
    "ncea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.5deg__UHAM-ICDC.nc')\n",
    "nc2ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_Bor_MAR_cru_ts3.22.1998.2010.pre.dat.nc')\n",
    "nc3ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_seasonal_soilw.mon.mean.1998-2010.v2.nc')\n",
    "#nc3ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_shift_min_bor_yearmean_soilw.mon.mean.1998-2010.v2.nc')\n",
    "nc4ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_0_bor_M6H_0.5_tmin.2m.gauss.1998-2010.nc')\n",
    "nc5ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_bor_GLC2000_0.5.nc')\n",
    "nc6ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_0_bor_0.5_PZI.nc')\n",
    "nc7ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_masked_bor_0.5_Year_Fire_Frequency_1996-2012.nc')\n",
    "# nc8ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_shift_bor_0.5_Burned_1996-2012_SUM.nc')\n",
    "nc9ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_shift_Bor_ReBil_GloElev_05.nc')\n",
    "nc10ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_Boreal_Mean_GDD0_1998-2010_0.5.nc')\n",
    "nc11ea_east = netCDF4.Dataset('/path_to_dataset/NEA_East_50_Bor_Ease25_Northen_Emisphere_Mean_Thawing_Depth_1901-2002.nc')\n",
    "nc12ea_east = netCDF4.Dataset('/Users/Beniamino/Documents/Datasets/Soil_texture/NEA_East_50_soil_texture_0.5.nc')\n",
    "nc13ea_east = netCDF4.Dataset('/path_to_dataset/Terra_NEA_East_50_MODIS_0.5_C6_GSSNDVI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "#nc13ea_east = netCDF4.Dataset('/path_to_dataset/Terra_NEA_East_50_MODIS_0.5_C6_GSSLAI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "\n",
    "#nc005ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.05deg__UHAM-ICDC.nc')\n",
    "ncea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.5deg__UHAM-ICDC.nc')\n",
    "nc2ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_Bor_MAR_cru_ts3.22.1998.2010.pre.dat.nc')\n",
    "nc3ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_seasonal_soilw.mon.mean.1998-2010.v2.nc')\n",
    "#nc3ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_shift_min_bor_yearmean_soilw.mon.mean.1998-2010.v2.nc')\n",
    "nc4ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_0_bor_M6H_0.5_tmin.2m.gauss.1998-2010.nc')\n",
    "nc5ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_bor_GLC2000_0.5.nc')\n",
    "nc6ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_0_bor_0.5_PZI.nc')\n",
    "nc7ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_masked_bor_0.5_Year_Fire_Frequency_1996-2012.nc')\n",
    "# nc8ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_shift_bor_0.5_Burned_1996-2012_SUM.nc')\n",
    "nc9ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_shift_Bor_ReBil_GloElev_05.nc')\n",
    "nc10ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_Boreal_Mean_GDD0_1998-2010_0.5.nc')\n",
    "nc11ea_west = netCDF4.Dataset('/path_to_dataset/NEA_West_Russia_50_Bor_Ease25_Northen_Emisphere_Mean_Thawing_Depth_1901-2002.nc')\n",
    "nc12ea_west = netCDF4.Dataset('/Users/Beniamino/Documents/Datasets/Soil_texture/NEA_West_Russia_50_soil_texture_0.5.nc')\n",
    "nc13ea_west = netCDF4.Dataset('/path_to_dataset/Terra_NEA_West_Russia_50_MODIS_0.5_C6_GSSNDVI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "#nc13ea_west = netCDF4.Dataset('/path_to_dataset/Terra_NEA_West_Russia_50_MODIS_0.5_C6_GSSLAI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "\n",
    "lon_ea_east = ncea_east.variables['lon']#.tolist()\n",
    "lon_ea_west = ncea_west.variables['lon']\n",
    "\n",
    "lat_ea_east = ncea_east.variables['lat']\n",
    "lat_ea_west = ncea_west.variables['lat']\n",
    "\n",
    "\n",
    "#lon005ea_east = nc005ea_east.variables['lon']#.tolist()\n",
    "#lon005ea_west = nc005ea_west.variables['lon']\n",
    "\n",
    "#lat005ea_east = nc005ea_east.variables['lat']\n",
    "\n",
    "\n",
    "vcff_ea_east = ncea_east.variables['vcff'][0]\n",
    "#vcff005_ea_east = nc005ea_east.variables['vcff'][0]\n",
    "pcp_ea_east = nc2ea_east.variables['pre'][0]\n",
    "sm_ea_east = nc3ea_east.variables['soilw']\n",
    "tm_ea_east = nc4ea_east.variables['tmin'][0]-273.15\n",
    "glc_ea_east = nc5ea_east.variables['var'][:]\n",
    "pzi_ea_east = nc6ea_east.variables['pzi'][:]\n",
    "ff_ea_east = nc7ea_east.variables['firefreq'][:]\n",
    "# ba_ea_east = nc8ea_east.variables['burnedarea_east'][:]\n",
    "elev_ea_east = nc9ea_east.variables['elev'][:]\n",
    "gdd0_ea_east = nc10ea_east.variables['gdd0'][0]\n",
    "ptd_ea_east = nc11ea_east.variables['thw_depth'][0]\n",
    "st_ea_east = nc12ea_east.variables['var600']\n",
    "lai_ea_east = nc13ea_east.variables['t_trend'][:]\n",
    "\n",
    "vcff_ea_west = ncea_west.variables['vcff'][0]\n",
    "#vcff005_ea_west = nc005ea_west.variables['vcff'][0]\n",
    "pcp_ea_west = nc2ea_west.variables['pre'][0]\n",
    "sm_ea_west = nc3ea_west.variables['soilw']\n",
    "tm_ea_west = nc4ea_west.variables['tmin'][0]-273.15\n",
    "glc_ea_west = nc5ea_west.variables['var'][:]\n",
    "pzi_ea_west = nc6ea_west.variables['pzi'][:]\n",
    "ff_ea_west = nc7ea_west.variables['firefreq'][:]\n",
    "# ba_ea_west = nc8ea_west.variables['burnedarea'][:]\n",
    "elev_ea_west = nc9ea_west.variables['elev'][:]\n",
    "gdd0_ea_west = nc10ea_west.variables['gdd0'][0]\n",
    "ptd_ea_west = nc11ea_west.variables['thw_depth'][0]\n",
    "st_ea_west = nc12ea_west.variables['var600']\n",
    "lai_ea_west = nc13ea_west.variables['t_trend'][:]\n",
    "\n",
    "\n",
    "\n",
    "picctr = 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create new variable of state accordind to the results of the optimisation\n",
    "vcf_state_ea_east = np.empty(shape=vcff_ea_east.shape)\n",
    "for i in range(vcff_ea_east.shape[0]):\n",
    "    for j in range(vcff_ea_east.shape[1]):\n",
    "        if vcff_ea_east[i][j]<0: \n",
    "            vcf_state_ea_east[i][j] = 0  # null\n",
    "        elif vcff_ea_east[i][j]<10:\n",
    "            vcf_state_ea_east[i][j] = 1  # bare, the first peak\n",
    "        elif vcff_ea_east[i][j]<20:\n",
    "            vcf_state_ea_east[i][j] = 1  # bare\n",
    "        elif vcff_ea_east[i][j]<45:\n",
    "            vcf_state_ea_east[i][j] = 3  # savanna, the third peak\n",
    "        elif vcff_ea_east[i][j]<100:\n",
    "            vcf_state_ea_east[i][j] = 4  # forest, the fourth peak\n",
    "        else:\n",
    "            vcf_state_ea_east[i][j] = 5  # null\n",
    "\n",
    "\n",
    "            # create new variable of state accordind to the results of the optimisation\n",
    "vcf_state_ea_west = np.empty(shape=vcff_ea_west.shape)\n",
    "for i in range(vcff_ea_west.shape[0]):\n",
    "    for j in range(vcff_ea_west.shape[1]):\n",
    "        if vcff_ea_west[i][j]<0: \n",
    "            vcf_state_ea_west[i][j] = 0  # null\n",
    "        elif vcff_ea_west[i][j]<10:\n",
    "            vcf_state_ea_west[i][j] = 1  # bare, the first peak\n",
    "        elif vcff_ea_west[i][j]<20:\n",
    "            vcf_state_ea_west[i][j] = 1  # bare\n",
    "        elif vcff_ea_west[i][j]<45:\n",
    "            vcf_state_ea_west[i][j] = 3  # savanna, the third peak\n",
    "        elif vcff_ea_west[i][j]<100:\n",
    "            vcf_state_ea_west[i][j] = 4  # forest, the fourth peak\n",
    "        else:\n",
    "            vcf_state_ea_west[i][j] = 5  # null\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Let's create a dataset with raveled arrays of data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# no need to execute this egain\n",
    "# create dataset for R with the raveled arrays of data\n",
    "\n",
    "## IMPORTANT: changes in numpy -> do np.ma.getdata(np.ravel())\n",
    "\n",
    "lons_ea_east,lats_ea_east = np.ma.getdata(np.meshgrid(lon_ea_east,lat_ea_east))\n",
    "arlons_ea_east, arlats_ea_east = np.ma.getdata(np.ravel(lons_ea_east)), np.ma.getdata(np.ravel(lats_ea_east))\n",
    "arsm_ea_east = np.ma.getdata(np.ravel(sm_ea_east[1]))\n",
    "arvcff_ea_east = np.ma.getdata(np.ravel(vcff_ea_east))\n",
    "arpcp_ea_east = np.ma.getdata(np.ravel(pcp_ea_east))\n",
    "artm_ea_east = np.ma.getdata(np.ravel(tm_ea_east))\n",
    "arglc_ea_east = np.ma.getdata(np.ravel(glc_ea_east))\n",
    "arpzi_ea_east = np.ma.getdata(np.ravel(pzi_ea_east))\n",
    "arff_ea_east = np.ma.getdata(np.ravel(ff_ea_east))\n",
    "#arba_ea_east = np.ravel(ba_ea_east)\n",
    "arelev_ea_east = np.ma.getdata(np.ravel(elev_ea_east))\n",
    "argdd0_ea_east = np.ma.getdata(np.ravel(gdd0_ea_east))\n",
    "arptd_ea_east = np.ma.getdata(np.ravel(ptd_ea_east))\n",
    "arstate_ea_east = np.ma.getdata(np.ravel(vcf_state_ea_east))\n",
    "arst_ea_east = np.ma.getdata(np.ravel(st_ea_east))\n",
    "arlai_ea_east = np.ma.getdata(np.ravel(lai_ea_east))\n",
    "\n",
    "\n",
    "count = 0\n",
    "index = []\n",
    "for i in range(len(arsm_ea_east)): \n",
    "    if (arsm_ea_east[i]<0 or arpcp_ea_east[i]>4000 or arvcff_ea_east[i]>100 or arglc_ea_east[i]>14 \\\n",
    "        or arpzi_ea_east[i]<0 or arff_ea_east[i]<0 or arelev_ea_east[i]>2000 or arelev_ea_east[i]<0 or argdd0_ea_east[i]<0 or \\\n",
    "       arptd_ea_east[i]<0):\n",
    "        index.append(i)\n",
    "        count+=1\n",
    "data_sm_ea_east = np.delete(arsm_ea_east,index)\n",
    "data_vcff_ea_east = np.delete(arvcff_ea_east,index)\n",
    "data_pcp_ea_east = np.delete(arpcp_ea_east,index)\n",
    "data_tm_ea_east = np.delete(artm_ea_east, index)\n",
    "data_pzi_ea_east = np.delete(arpzi_ea_east, index)\n",
    "data_ff_ea_east = np.delete(arff_ea_east, index)\n",
    "#data_ba_ea_east = np.delete(arba, index)\n",
    "data_elev_ea_east = np.delete(arelev_ea_east, index)\n",
    "data_gdd0_ea_east = np.delete(argdd0_ea_east, index)\n",
    "data_ptd_ea_east = np.delete(arptd_ea_east, index)\n",
    "data_state_ea_east = np.delete(arstate_ea_east, index)\n",
    "data_st_ea_east = np.delete(arst_ea_east, index)\n",
    "data_lai_ea_east = np.delete(arlai_ea_east, index)\n",
    "\n",
    "\n",
    "data_lons_ea_east, data_lats_ea_east = np.delete(arlons_ea_east,index), np.delete(arlats_ea_east,index)\n",
    "\n",
    "\n",
    "#dataforr_ea_east = np.transpose(np.array([data_vcff_ea_east, data_pcp_ea_east, data_sm_ea_east, data_tm_ea_east, data_pzi_ea_east,\\\n",
    "#                                     data_ff_ea_east, data_gdd0_ea_east, data_ptd_ea_east, data_st_ea_east, data_state_ea_east, data_lai_ea_east, data_elev_ea_east, data_lons_ea_east, data_lats_ea_east])) \n",
    "#dataframe_ea_east_all = DataFrame(dataforr_ea_east[:,0:10], columns=[\"VCFF\", \"MAR\", \"MASM\", \"Tmin\", \"PZI\", \"FF\", \"GDD0\", \"PTD\", \"ST\",\"State\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# save the data to a txt that can be reloaded  \n",
    "# it is executed twice, once including lai, once including ndvi\n",
    "lai_ea_east = np.loadtxt('boreal_ea_east_database_lai_included.txt')[:,10]  #what is here called lai is actually ndvi, whereas this lai_na_east is the actual lai\n",
    "dataforr_ea_east = np.transpose(np.array([data_vcff_ea_east, data_pcp_ea_east, data_sm_ea_east, data_tm_ea_east, data_pzi_ea_east,\\\n",
    "                                     data_ff_ea_east, data_gdd0_ea_east, data_ptd_ea_east, data_st_ea_east, data_state_ea_east, lai_ea_east, data_lai_ea_east, data_elev_ea_east, data_lons_ea_east, data_lats_ea_east])) \n",
    "\n",
    "np.savetxt(\"boreal_ea_east_database_lai_ndvi_included.txt\", dataforr_ea_east)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# no need to execute this egain\n",
    "# create dataset for R with the raveled arrays of data\n",
    "\n",
    "## IMPORTANT: changes in numpy -> do np.ma.getdata(np.ravel())\n",
    "\n",
    "lons_ea_west,lats_ea_west = np.ma.getdata(np.meshgrid(lon_ea_west,lat_ea_west))\n",
    "arlons_ea_west, arlats_ea_west = np.ma.getdata(np.ravel(lons_ea_west)), np.ma.getdata(np.ravel(lats_ea_west))\n",
    "arsm_ea_west = np.ma.getdata(np.ravel(sm_ea_west[1]))\n",
    "arvcff_ea_west = np.ma.getdata(np.ravel(vcff_ea_west))\n",
    "arpcp_ea_west = np.ma.getdata(np.ravel(pcp_ea_west))\n",
    "artm_ea_west = np.ma.getdata(np.ravel(tm_ea_west))\n",
    "arglc_ea_west = np.ma.getdata(np.ravel(glc_ea_west))\n",
    "arpzi_ea_west = np.ma.getdata(np.ravel(pzi_ea_west))\n",
    "arff_ea_west = np.ma.getdata(np.ravel(ff_ea_west))\n",
    "#arba_ea_west = np.ravel(ba_ea_west)\n",
    "arelev_ea_west = np.ma.getdata(np.ravel(elev_ea_west))\n",
    "argdd0_ea_west = np.ma.getdata(np.ravel(gdd0_ea_west))\n",
    "arptd_ea_west = np.ma.getdata(np.ravel(ptd_ea_west))\n",
    "arstate_ea_west = np.ma.getdata(np.ravel(vcf_state_ea_west))\n",
    "arst_ea_west = np.ma.getdata(np.ravel(st_ea_west))\n",
    "arlai_ea_west = np.ma.getdata(np.ravel(lai_ea_west))\n",
    "\n",
    "\n",
    "count = 0\n",
    "index = []\n",
    "for i in range(len(arsm_ea_west)): \n",
    "    if (arsm_ea_west[i]<0 or arpcp_ea_west[i]>4000 or arvcff_ea_west[i]>100 or arglc_ea_west[i]>14 \\\n",
    "        or arpzi_ea_west[i]<0 or arff_ea_west[i]<0 or arelev_ea_west[i]>2000 or arelev_ea_west[i]<0 or argdd0_ea_west[i]<0 or \\\n",
    "       arptd_ea_west[i]<0):\n",
    "        index.append(i)\n",
    "        count+=1\n",
    "data_sm_ea_west = np.delete(arsm_ea_west,index)\n",
    "data_vcff_ea_west = np.delete(arvcff_ea_west,index)\n",
    "data_pcp_ea_west = np.delete(arpcp_ea_west,index)\n",
    "data_tm_ea_west = np.delete(artm_ea_west, index)\n",
    "data_pzi_ea_west = np.delete(arpzi_ea_west, index)\n",
    "data_ff_ea_west = np.delete(arff_ea_west, index)\n",
    "#data_ba_ea_west = np.delete(arba, index)\n",
    "data_elev_ea_west = np.delete(arelev_ea_west, index)\n",
    "data_gdd0_ea_west = np.delete(argdd0_ea_west, index)\n",
    "data_ptd_ea_west = np.delete(arptd_ea_west, index)\n",
    "data_state_ea_west = np.delete(arstate_ea_west, index)\n",
    "data_st_ea_west = np.delete(arst_ea_west, index)\n",
    "data_lai_ea_west = np.delete(arlai_ea_west, index)\n",
    "\n",
    "data_lons_ea_west, data_lats_ea_west = np.delete(arlons_ea_west,index), np.delete(arlats_ea_west,index)\n",
    "\n",
    "\n",
    "#dataframe_ea_west_all = DataFrame(dataforr_ea_west[:,0:9], columns=[\"VCFF\", \"MAR\", \"MASM\", \"Tmin\", \"PZI\", \"FF\", \"GDD0\", \"PTD\",\"State\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# save the data to a txt that can be reloaded  \n",
    "\n",
    "lai_ea_west = np.loadtxt('boreal_ea_west_database_lai_included.txt')[:,10]  #what is here called lai is actually ndvi, whereas this lai_na_east is the actual lai\n",
    "\n",
    "dataforr_ea_west = np.transpose(np.array([data_vcff_ea_west, data_pcp_ea_west, data_sm_ea_west, data_tm_ea_west, data_pzi_ea_west,\\\n",
    "                                     data_ff_ea_west, data_gdd0_ea_west, data_ptd_ea_west, data_st_ea_west, data_state_ea_west, lai_ea_west, data_lai_ea_west, data_elev_ea_west, data_lons_ea_west, data_lats_ea_west])) \n",
    "\n",
    "np.savetxt(\"boreal_ea_west_database_lai_ndvi_included.txt\", dataforr_ea_west)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# from here on it is actually north america\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# execute this only if you want to do the basemap plot\n",
    "# pah_to_dataset = Users/Beniamino/Documents/Datasets/Boreal\n",
    "nc005ea = netCDF4.Dataset('/path_to_dataset/\\\n",
    "NA_West_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.05deg__UHAM-ICDC.nc')\n",
    "ncea = netCDF4.Dataset('/path_to_dataset/\\\n",
    "NA_West_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.5deg__UHAM-ICDC.nc')\n",
    "nc2ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_Bor_MAR_cru_ts3.22.1998.2010.pre.dat.nc')\n",
    "nc3ea = netCDF4.Dataset('/path_to_dataset/NA_West_seasonal_soilw.mon.mean.1998-2010.v2.nc')\n",
    "#nc3ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_min_bor_yearmean_soilw.mon.mean.1998-2010.v2.nc')\n",
    "nc4ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_0_bor_M6H_0.5_tmin.2m.gauss.1998-2010.nc')\n",
    "nc5ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_bor_GLC2000_0.5.nc')\n",
    "nc6ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_0_bor_0.5_PZI.nc')\n",
    "nc7ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_masked_bor_0.5_Year_Fire_Frequency_1996-2012.nc')\n",
    "# nc8ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_bor_0.5_Burned_1996-2012_SUM.nc')\n",
    "nc9ea = netCDF4.Dataset('/path_to_dataset/NA_West_shift_Bor_ReBil_GloElev_05.nc')\n",
    "nc10ea = netCDF4.Dataset('/path_to_dataset/NA_West_Boreal_Mean_GDD0_1998-2010_0.5.nc')\n",
    "nc11ea = netCDF4.Dataset('/path_to_dataset/NA_West_Bor_Ease25_Northen_Emisphere_Mean_Thawing_Depth_1901-2002.nc')\n",
    "nc12ea = netCDF4.Dataset('/Users/Beniamino/Documents/Datasets/Soil_texture/NA_West_soil_texture_0.5.nc')\n",
    "#nc13ea = netCDF4.Dataset('/path_to_dataset/Terra_NA_West_MODIS_0.5_C6_GSSLAI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "nc13ea = netCDF4.Dataset('/path_to_dataset/Terra_NA_West_MODIS_0.5_C6_GSSNDVI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "\n",
    "nc005na = netCDF4.Dataset('/path_to_dataset/\\\n",
    "NA_East_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.05deg__UHAM-ICDC.nc')\n",
    "ncna = netCDF4.Dataset('/path_to_dataset/\\\n",
    "NA_East_shift_Bor_MOD44B_V5__VCFFplusVCF__2010__0.5deg__UHAM-ICDC.nc')\n",
    "nc2na = netCDF4.Dataset('/path_to_dataset/NA_East_shift_Bor_MAR_cru_ts3.22.1998.2010.pre.dat.nc')\n",
    "nc3na = netCDF4.Dataset('/path_to_dataset/NA_East_seasonal_soilw.mon.mean.1998-2010.v2.nc')\n",
    "#nc3na = netCDF4.Dataset('/path_to_dataset/NA_shift_min_bor_yearmean_soilw.mon.mean.1998-2010.v2.nc')\n",
    "nc4na = netCDF4.Dataset('/path_to_dataset/NA_East_shift_0_bor_M6H_0.5_tmin.2m.gauss.1998-2010.nc')\n",
    "nc5na = netCDF4.Dataset('/path_to_dataset/NA_East_shift_bor_GLC2000_0.5.nc')\n",
    "nc6na = netCDF4.Dataset('/path_to_dataset/NA_East_shift_0_bor_0.5_PZI.nc')\n",
    "nc7na = netCDF4.Dataset('/path_to_dataset/NA_East_shift_masked_bor_0.5_Year_Fire_Frequency_1996-2012.nc')\n",
    "# nc8na = netCDF4.Dataset('/path_to_dataset/NA_shift_bor_0.5_Burned_1996-2012_SUM.nc')\n",
    "nc9na = netCDF4.Dataset('/path_to_dataset/NA_East_shift_Bor_ReBil_GloElev_05.nc')\n",
    "nc10na = netCDF4.Dataset('/path_to_dataset/NA_East_Boreal_Mean_GDD0_1998-2010_0.5.nc')\n",
    "nc11na = netCDF4.Dataset('/path_to_dataset/NA_East_Bor_Ease25_Northen_Emisphere_Mean_Thawing_Depth_1901-2002.nc')\n",
    "nc12na = netCDF4.Dataset('/Users/Beniamino/Documents/Datasets/Soil_texture/NA_East_soil_texture_0.5.nc')\n",
    "#nc13na = netCDF4.Dataset('/path_to_dataset/Terra_NA_East_MODIS_0.5_C6_GSSLAI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "nc13na = netCDF4.Dataset('/path_to_dataset/Terra_NA_East_MODIS_0.5_C6_GSSNDVI_2000_2015_Decadal_Trend_SigOnly.nc')\n",
    "\n",
    "\n",
    "lon_ea = ncea.variables['lon']#.tolist()\n",
    "lon_na = ncna.variables['lon']\n",
    "\n",
    "lat = ncea.variables['lat']\n",
    "\n",
    "lon005ea = nc005ea.variables['lon']#.tolist()\n",
    "lon005na = nc005na.variables['lon']\n",
    "\n",
    "lat005ea = nc005ea.variables['lat']\n",
    "\n",
    "\n",
    "vcff_ea = ncea.variables['vcff'][0]\n",
    "vcff005_ea = nc005ea.variables['vcff'][0]\n",
    "\n",
    "pcp_ea = nc2ea.variables['pre'][0]\n",
    "sm_ea = nc3ea.variables['soilw']\n",
    "tm_ea = nc4ea.variables['tmin'][0]-273.15\n",
    "glc_ea = nc5ea.variables['var'][:]\n",
    "pzi_ea = nc6ea.variables['pzi'][:]\n",
    "ff_ea = nc7ea.variables['firefreq'][:]\n",
    "# ba_ea = nc8ea.variables['burnedarea'][:]\n",
    "elev_ea = nc9ea.variables['elev'][:]\n",
    "gdd0_ea = nc10ea.variables['gdd0'][0]\n",
    "st_ea = nc12ea.variables['var600']\n",
    "ptd_ea = nc11ea.variables['thw_depth'][0]\n",
    "lai_ea = nc13ea.variables['t_trend'][:]\n",
    "\n",
    "\n",
    "vcff_na = ncna.variables['vcff'][0]\n",
    "vcff005_na = nc005na.variables['vcff'][0]\n",
    "\n",
    "pcp_na = nc2na.variables['pre'][0]\n",
    "sm_na = nc3na.variables['soilw']\n",
    "tm_na = nc4na.variables['tmin'][0]-273.15\n",
    "glc_na = nc5na.variables['var'][:]\n",
    "pzi_na = nc6na.variables['pzi'][:]\n",
    "ff_na = nc7na.variables['firefreq'][:]\n",
    "# ba_na = nc8na.variables['burnedarea'][:]\n",
    "elev_na = nc9na.variables['elev'][:]\n",
    "gdd0_na = nc10na.variables['gdd0'][0]\n",
    "ptd_na = nc11na.variables['thw_depth'][0]\n",
    "lai_na = nc13na.variables['t_trend'][:]\n",
    "st_na = nc12na.variables['var600']\n",
    "\n",
    "\n",
    "picctr = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# create new variable of state accordind to the results of the optimisation\n",
    "vcf_state_na = np.empty(shape=vcff_na.shape)\n",
    "for i in range(vcff_na.shape[0]):\n",
    "    for j in range(vcff_na.shape[1]):\n",
    "        if vcff_na[i][j]<0: \n",
    "            vcf_state_na[i][j] = 0  # null\n",
    "        elif vcff_na[i][j]<10:\n",
    "            vcf_state_na[i][j] = 1  # bare, the first peak\n",
    "        elif vcff_na[i][j]<20:\n",
    "            vcf_state_na[i][j] = 1  # sparse savanna, the second peak\n",
    "        elif vcff_na[i][j]<45:\n",
    "            vcf_state_na[i][j] = 3  # savanna, the third peak\n",
    "        elif vcff_na[i][j]<100:\n",
    "            vcf_state_na[i][j] = 4  # forest, the fourth peak\n",
    "        else:\n",
    "            vcf_state_na[i][j] = 5  # null\n",
    "            \n",
    "            \n",
    "# create new variable of state accordind to the results of the optimisation\n",
    "vcf_state_ea = np.empty(shape=vcff_ea.shape)\n",
    "for i in range(vcff_ea.shape[0]):\n",
    "    for j in range(vcff_ea.shape[1]):\n",
    "        if vcff_ea[i][j]<0: \n",
    "            vcf_state_ea[i][j] = 0  # null\n",
    "        elif vcff_ea[i][j]<10:\n",
    "            vcf_state_ea[i][j] = 1  # bare, the first peak\n",
    "        elif vcff_ea[i][j]<20:\n",
    "            vcf_state_ea[i][j] = 1  # sparse savanna, the second peak\n",
    "        elif vcff_ea[i][j]<45:\n",
    "            vcf_state_ea[i][j] = 3  # savanna, the third peak\n",
    "        elif vcff_ea[i][j]<100:\n",
    "            vcf_state_ea[i][j] = 4  # forest, the fourth peak\n",
    "        else:\n",
    "            vcf_state_ea[i][j] = 5  # null\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# no need to execute this egain\n",
    "# create dataset for R with the raveled arrays of data\n",
    "\n",
    "## IMPORTANT: changes in numpy -> do np.ma.getdata(np.ravel())\n",
    "\n",
    "lons,lats = np.ma.getdata(np.meshgrid(lon_na,lat))\n",
    "arlons, arlats = np.ma.getdata(np.ravel(lons)), np.ma.getdata(np.ravel(lats))\n",
    "arsm_na = np.ma.getdata(np.ravel(sm_na[1]))\n",
    "arvcff_na = np.ma.getdata(np.ravel(vcff_na))\n",
    "arpcp_na = np.ma.getdata(np.ravel(pcp_na))\n",
    "artm_na = np.ma.getdata(np.ravel(tm_na))\n",
    "arglc_na = np.ma.getdata(np.ravel(glc_na))\n",
    "arpzi_na = np.ma.getdata(np.ravel(pzi_na))\n",
    "arff_na = np.ma.getdata(np.ravel(ff_na))\n",
    "#arba_na = np.ravel(ba_na)\n",
    "arelev_na = np.ma.getdata(np.ravel(elev_na))\n",
    "argdd0_na = np.ma.getdata(np.ravel(gdd0_na))\n",
    "arptd_na = np.ma.getdata(np.ravel(ptd_na))\n",
    "arstate_na = np.ma.getdata(np.ravel(vcf_state_na))\n",
    "arst_na = np.ma.getdata(np.ravel(st_na))\n",
    "arlai_na = np.ma.getdata(np.ravel(lai_na))\n",
    "\n",
    "\n",
    "count = 0\n",
    "index = []\n",
    "for i in range(len(arsm_na)): \n",
    "    if (arsm_na[i]<0 or arpcp_na[i]>4000 or arvcff_na[i]>100 or arglc_na[i]>14 \\\n",
    "        or arpzi_na[i]<0 or arff_na[i]<0 or arelev_na[i]>2000 or arelev_na[i]<0 or argdd0_na[i]<0 or \\\n",
    "       arptd_na[i]<0):\n",
    "        index.append(i)\n",
    "        count+=1\n",
    "data_sm_na = np.delete(arsm_na,index)\n",
    "data_vcff_na = np.delete(arvcff_na,index)\n",
    "data_pcp_na = np.delete(arpcp_na,index)\n",
    "data_tm_na = np.delete(artm_na, index)\n",
    "data_pzi_na = np.delete(arpzi_na, index)\n",
    "data_ff_na = np.delete(arff_na, index)\n",
    "#data_ba_na = np.delete(arba, index)\n",
    "data_elev_na = np.delete(arelev_na, index)\n",
    "data_gdd0_na = np.delete(argdd0_na, index)\n",
    "data_ptd_na = np.delete(arptd_na, index)\n",
    "data_state_na = np.delete(arstate_na, index)\n",
    "data_st_na = np.delete(st_na, index)\n",
    "data_lai_na = np.delete(arlai_na, index)\n",
    "\n",
    "\n",
    "\n",
    "data_lons, data_lats = np.delete(arlons,index), np.delete(arlats,index)\n",
    "\n",
    "\n",
    "#dataforr_na = np.transpose(np.array([data_vcff_na, data_pcp_na, data_sm_na, data_tm_na, data_pzi_na,\\\n",
    "#                                     data_ff_na, data_gdd0_na, data_ptd_na, data_st_na, data_state_na, data_lai_na, data_elev_na, data_lons, data_lats])) \n",
    "dataframe_na_all = DataFrame(dataforr_na[:,0:10], columns=[\"VCFF\", \"MAR\", \"MASM\", \"Tmin\", \"PZI\", \"FF\", \"GDD0\", \"PTD\", \"ST\",\"State\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "canadianupdate_na_east=np.loadtxt(\"boreal_na_east_database_canadianupdate_without_soiltexture.txt\")\n",
    "lai_na_east = np.loadtxt('boreal_na_east_database_lai_included.txt')[:,10]  #what is here called lai is actually ndvi, whereas this lai_na_east is the actual lai\n",
    "dataforr_na = np.transpose(np.array([data_vcff_na, data_pcp_na, data_sm_na, data_tm_na, data_pzi_na,canadianupdate_na_east[:,5], \\\n",
    "                                     data_gdd0_na, data_ptd_na, data_st_na, data_state_na, lai_na_east, data_lai_na, data_elev_na, data_lons, data_lats])) \n",
    "np.savetxt(\"boreal_na_east_database_lai_ndvi_included.txt\", dataforr_na)\n",
    "#np.savetxt(\"boreal_na_east_database_lai_included.txt\", dataforr_na)\n",
    "#np.savetxt(\"boreal_na_east_database_canadianupdate.txt\", dataforr_na)\n",
    "# save the data to a txt that can be reloaded  \n",
    "#np.savetxt(\"boreal_na_east_database_canadianupdate_without_soiltexture.txt\", canadianupdate_na_east)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# no need to execute this egain\n",
    "# create dataset for R with the raveled arrays of data\n",
    "\n",
    "## IMPORTANT: changes in numpy -> do np.ma.getdata(np.ravel())\n",
    "\n",
    "lons,lats = np.ma.getdata(np.meshgrid(lon_ea,lat))\n",
    "arlons, arlats = np.ma.getdata(np.ravel(lons)), np.ma.getdata(np.ravel(lats))\n",
    "arsm_ea = np.ma.getdata(np.ravel(sm_ea[1]))\n",
    "arvcff_ea = np.ma.getdata(np.ravel(vcff_ea))\n",
    "arpcp_ea = np.ma.getdata(np.ravel(pcp_ea))\n",
    "artm_ea = np.ma.getdata(np.ravel(tm_ea))\n",
    "arglc_ea = np.ma.getdata(np.ravel(glc_ea))\n",
    "arpzi_ea = np.ma.getdata(np.ravel(pzi_ea))\n",
    "arff_ea = np.ma.getdata(np.ravel(ff_ea))\n",
    "#arba_ea = np.ravel(ba_ea)\n",
    "arelev_ea = np.ma.getdata(np.ravel(elev_ea))\n",
    "argdd0_ea = np.ma.getdata(np.ravel(gdd0_ea))\n",
    "arptd_ea = np.ma.getdata(np.ravel(ptd_ea))\n",
    "arstate_ea = np.ma.getdata(np.ravel(vcf_state_ea))\n",
    "arst_ea = np.ma.getdata(np.ravel(st_ea))\n",
    "arlai_ea = np.ma.getdata(np.ravel(lai_ea))\n",
    "\n",
    "\n",
    "count = 0\n",
    "index = []\n",
    "for i in range(len(arsm_ea)): \n",
    "    if (arsm_ea[i]<0 or arpcp_ea[i]>4000 or arvcff_ea[i]>100 or arglc_ea[i]>14 \\\n",
    "        or arpzi_ea[i]<0 or arff_ea[i]<0 or arelev_ea[i]>2000 or arelev_ea[i]<0 or argdd0_ea[i]<0 or \\\n",
    "       arptd_ea[i]<0):\n",
    "        index.append(i)\n",
    "        count+=1\n",
    "data_sm_ea = np.delete(arsm_ea,index)\n",
    "data_vcff_ea = np.delete(arvcff_ea,index)\n",
    "data_pcp_ea = np.delete(arpcp_ea,index)\n",
    "data_tm_ea = np.delete(artm_ea, index)\n",
    "data_pzi_ea = np.delete(arpzi_ea, index)\n",
    "data_ff_ea = np.delete(arff_ea, index)\n",
    "#data_ba_ea = np.delete(arba, index)\n",
    "data_elev_ea = np.delete(arelev_ea, index)\n",
    "data_gdd0_ea = np.delete(argdd0_ea, index)\n",
    "data_ptd_ea = np.delete(arptd_ea, index)\n",
    "data_state_ea = np.delete(arstate_ea, index)\n",
    "data_st_ea = np.delete(st_ea, index)\n",
    "data_lai_ea = np.delete(arlai_ea, index)\n",
    "\n",
    "\n",
    "data_lons, data_lats = np.delete(arlons,index), np.delete(arlats,index)\n",
    "\n",
    "\n",
    "#dataforr_ea = np.transpose(np.array([data_vcff_ea, data_pcp_ea, data_sm_ea, data_tm_ea, data_pzi_ea,\\\n",
    "#                                     data_ff_ea, data_gdd0_ea, data_ptd_ea, data_st_ea, data_state_ea, data_lai_ea, data_elev_ea, data_lons, data_lats])) \n",
    "#dataframe_ea_all = DataFrame(dataforr_ea[:,0:10], columns=[\"VCFF\", \"MAR\", \"MASM\", \"Tmin\", \"PZI\", \"FF\", \"GDD0\", \"PTD\", \"ST\", \"State\"])\n",
    "\n",
    "# save the data to a txt that can be reloaded  \n",
    "#np.savetxt(\"boreal_na_west_database.txt\", dataforr_ea)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "canadianupdate_na_west=np.loadtxt(\"boreal_na_west_database_canadianupdate_without_soiltexture.txt\")\n",
    "lai_na_west = np.loadtxt('boreal_na_west_database_lai_included.txt')[:,10]  #what is here called lai is actually ndvi, whereas this lai_na_east is the actual lai\n",
    "dataforr_ea = np.transpose(np.array([data_vcff_ea, data_pcp_ea, data_sm_ea, data_tm_ea, data_pzi_ea,\\\n",
    "                                     canadianupdate_na_west[:,5], data_gdd0_ea, data_ptd_ea, data_st_ea, data_state_ea, lai_na_west, data_lai_ea, data_elev_ea, data_lons, data_lats])) \n",
    "np.savetxt(\"boreal_na_west_database_lai_ndvi_included.txt\", dataforr_ea)\n",
    "#np.savetxt(\"boreal_na_west_database_canadianupdate.txt\", dataforr_ea)\n",
    "#np.savetxt(\"boreal_na_west_database_canadianupdate_without_soiltexture.txt\", canadianupdate_na_west)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}
