Large_Melt_Map.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Large_Melt_Map\n",
    "\n",
    "Calculate SIR and GRD melt onset dates by year for selected subset areas.\n",
    "\n",
    "Also calculates average MOD across years at each pixel.\n",
    "\n",
    "Saves SIR or GRD MOD data together with geolocation information in pickle files\n",
    "\n",
    "Makes maps of various annual and/or average results.\n",
    "\n",
    "FIXME: \n",
    "The pixel overlays are simply a scatter plot with rectangular markers and are not using EASE-grid projection\n",
    "\n",
    "FIXME:\n",
    "Add a location for where to find the RGI overlays\n",
    "and read/display the RGI shapefiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in all the modules needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# check if a windows machine, it needs special attention\n",
    "# this extra step will bypass an error from mpl_toolkits.basemap\n",
    "import os\n",
    "if os.name == 'nt':\n",
    "    os.environ[\"PROJ_LIB\"] = os.path.join(os.environ[\"CONDA_PREFIX\"], \"Library\", \"share\")\n",
    "    os.environ[\"GDAL_DATA\"] = os.path.join(os.environ[\"CONDA_PREFIX\"], \"Library\", \"share\", \"gdal\")\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from netCDF4 import Dataset, num2date\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import re\n",
    "from cetbtools.ease2conv import Ease2Transform\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import glob\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_rows', None)\n",
    "pd.set_option('display.min_rows', None)\n",
    "pd.set_option('display.max_columns', 20)\n",
    "pd.set_option('display.width', 120)\n",
    "pd.set_option('display.max_colwidth', 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the local machine location of CETB data cubes\n",
    "# This directory is expected to contain subdirectories in the following hierarchy\n",
    "# that duplicates the hierarchy on the Google Shared Drive NSIDC-SD-CETB/v1/, \n",
    "# for example:\n",
    "# dataDir/F13_SSMI/N/nc_cubes/cubes_<regionName>\n",
    "#\n",
    "# outDir is the location where pkl files with the MOD data frames and geotiffs\n",
    "# will be stored. We are keeping it separate from the large cubefile data store\n",
    "# so that individual people can keep track of their own output files separately\n",
    "user = 'Joan' #Mariah #MJWindows #MAHMac\n",
    "if ('Joan' == user):\n",
    "    #dataDir = '/mnt/data3/cetb/nsidc0630_v1/' #jmr machine fringe \n",
    "    dataDir = Path(Path.home(), 'ceph', 'jmr204group','CETB_cubes')\n",
    "    scriptDir = Path(Path.home(), 'ipynb_melt_onset', 'scripts')\n",
    "    pklDir = Path(dataDir, 'pkls')\n",
    "#    outDir = Path(Path.home(), 'cetb/ipynb_melt__onset_plots')  #may need a spot for derived data\n",
    "#    outDir = Path(Path.home(), 'rdrive', 'jmr204group', 'ramage', 'geotiffs')\n",
    "    outDir = Path(Path.home(), 'noisy')\n",
    "elif ('Mariah' == user):\n",
    "    # dataDir = Path('R:\\\\jmr204group\\CETB_cubes')  # Mariah's PC\n",
    "    # changing this syntax to fix warning about invalid character '\\C'\n",
    "    dataDir = Path('R:', 'jmr204group', 'CETB_cubes')\n",
    "    scriptDir = Path(Path.home(), 'ipynb_melt_onset', 'scripts')\n",
    "    outDir = Path(Path.home(), 'ipynb_melt_onset_plots') # may want to update this to a geotiff output directory\n",
    "elif ('MJWindows' == user):\n",
    "    dataDir = Path('Z:/mj On My Mac/nsidc0630_v1') # Mary Jo's Windows machine\n",
    "    scriptDir = Path(Path.home(), 'ipynb_melt_onset', 'scripts')\n",
    "    pklDir = Path(dataDir, 'pkls')\n",
    "    outDir = Path(Path.home(), 'ipynb_melt_onset_plots') # may want to update this to a geotiff output directory\n",
    "elif ('MJMac' == user):\n",
    "    dataDir = Path(Path.home(), 'nsidc0630_v1') # Mary Jo's Mac\n",
    "    scriptDir = Path(Path.home(), 'ipynb_melt_onset', 'scripts')  \n",
    "    pklDir = Path(Path.home(), 'nsidc0630_v1', 'pkls')\n",
    "    outDir = Path(Path.home(), 'nsidc0630_v1')\n",
    "elif ('MAHMac' == user):\n",
    "    dataDir = Path(Path.home(), 'nsidc0630_v1') # Mary Jo's Mac\n",
    "    scriptDir = Path(Path.home(), 'Projects', 'ipynb_melt_onset', 'scripts')  \n",
    "    pklDir = Path(Path.home(), 'nsidc0630_v1', 'pkls')\n",
    "    outDir = Path(Path.home(), 'nsidc0630_v1')\n",
    "else:\n",
    "    raise ValueError(\"unknown user= %s\\n\" % (user) )\n",
    "    \n",
    "%cd $scriptDir\n",
    "outDir, dataDir, user"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the custom functions\n",
    "from CETB_IO import read_Tb_whole\n",
    "from CETB_IO import coords\n",
    "from CETB_IO import get_site_boundaries\n",
    "from CETB_algorithms import calc_DAV\n",
    "from CETB_IO import find_cube_offset\n",
    "from CETB_IO import grid_locations_of_subset\n",
    "from CETB_IO import years_for\n",
    "from CETB_IO import get_sir_info\n",
    "from CETB_algorithms import DAV_MOD\n",
    "from CETB_analysis import MOD_array\n",
    "from CETB_algorithms import end_high_DAV\n",
    "from CETB_IO import read_Tb_std_dev\n",
    "from CETB_IO import read_Tb\n",
    "from CETB_IO import write_df_to_geotiff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify inputs\n",
    "\n",
    "This cell is used to specify region, satellite, sensor, channel, and image reconstruction algorithm\n",
    "values that will be used to find the filenames with the data we are intereted in\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region='AKYukon'  #'GLaIL'  #make this the same syntax as cubefilenames and sub-directory\n",
    "sat_GRD='F18'  #'AQUA' for AMSRE, 'F13','F14','F15'... for SSMI\n",
    "sat_SIR= 'F18'\n",
    "sensor_GRD='SSMIS'#'AMSRE', 'SSMI', etc.\n",
    "sensor_SIR='SSMIS'\n",
    "channel_GRD='37V'  #'36V','36H', '18V','18H', etc. '19V','19H' and '37V','37H' for SSMI)\n",
    "channel_SIR='37V'\n",
    "alg_GRD='GRD'   #SIR or GRD\n",
    "alg_SIR='SIR'\n",
    "\n",
    "hemName = 'N'    \n",
    "\n",
    "# get sir to grd factor and sir_gpd name\n",
    "sir_2_grd_factor, gpd_SIR = get_sir_info(channel_SIR, hem=hemName)\n",
    "gpd_GRD = 'EASE2_%s25km' % hemName\n",
    "\n",
    "cubeType_GRD = channel_GRD + '-' + alg_GRD\n",
    "cubeType_SIR = channel_SIR + '-' + alg_SIR\n",
    "\n",
    "# When we change to v2 files, there will be no need for different providers by sensor, since\n",
    "# there is no provider name in the v2 filenames\n",
    "if ('SSMI' == sensor_GRD) or ('SSMIS' == sensor_GRD):\n",
    "    provider='CSU' # FIXME, this will have to change for v2 files\n",
    "    version='v1.*'\n",
    "elif 'AMSRE' == sensor_GRD:\n",
    "    provider='RSS' # FIXME: this will have to change for v2 files--it \n",
    "    version='v1.3'\n",
    "\n",
    "\n",
    "# on Joan's machine\n",
    "#datadir_GRD = dataDir + sat_GRD+'_'+sensor_GRD+'/'+region+'/' \n",
    "#datadir_SIR = dataDir + sat_SIR+'_'+sensor_SIR+'/'+region+'/' \n",
    "# on MJ's machine\n",
    "datadir_GRD = \"%s/%s_%s/%s/nc_cubes/cubes_%s/\" % (\n",
    "    dataDir, sat_GRD, sensor_GRD, hemName, region )\n",
    "datadir_SIR = \"%s/%s_%s/%s/nc_cubes/cubes_%s/\" % (\n",
    "    dataDir, sat_SIR, sensor_SIR, hemName, region )\n",
    "\n",
    "# prefix filepath\n",
    "prefix_GRD = 'CETB.cubefile.'+region+'.'+sat_GRD+'_'+sensor_GRD+'-'+channel_GRD+'-'+alg_GRD+'-'+provider+'-'+version\n",
    "prefix_SIR = 'CETB.cubefile.'+region+'.'+sat_SIR+'_'+sensor_SIR+'-'+channel_SIR+'-'+alg_SIR+'-'+provider+'-'+version\n",
    "\n",
    "Years=years_for(sat_GRD)\n",
    "subYears=Years[0:2]\n",
    "#might want to truncate Years to subset if very slow during testing\n",
    "#if we give it more years than available what do we want it to do? \n",
    "#warn me but return what it finds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify geographic location that we are interested in processing\n",
    "\n",
    "Do lat and long in decimal degrees, choose corners that are completely outside the set of pixels that we want to process.\n",
    "'SiteLabel' should be short and not have spaces or punctuation.\n",
    "Set a 'Site' name for each area that will be used for plots/labels later on in this notebook. Site names can have spaces/punctuation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_site_boundaries('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SiteLabel='barrow'\n",
    "lat_start, lat_end, lon_start, lon_end, Site=get_site_boundaries(SiteLabel)\n",
    "lat_start, lat_end, lon_start, lon_end, Site"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get cube offset for finding row/col of the region's cube relative to the full EASE2 hemisphere grid\n",
    "# function is region specific\n",
    "# The returned value from find_cube_offset is the (row, col) in the EASE2_N 3.125 resolution grid that\n",
    "# corresponds to the center of the UL grid cell of the subset cube\n",
    "find_cube_offset(region, cubeDir=datadir_SIR, cubeType=cubeType_SIR, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the GRD pixel coordinates relative to the cube UL for the lat/lon rectangle chosen in our region\n",
    "# and then calculate the corrsponding SIR \"envelope\" pixel row/col numbers by multiplying by the grid factor relationship\n",
    "rows_cols_GRD = coords(datadir_GRD, prefix_GRD, lat_start, lat_end, lon_start, lon_end)\n",
    "rows_cols_env = tuple(np.array(rows_cols_GRD) * sir_2_grd_factor)\n",
    "print(rows_cols_GRD)\n",
    "print(rows_cols_env)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in SIR TB data\n",
    "# Truncate Years here for speed \n",
    "data_SIR = read_Tb_whole(datadir_SIR, prefix_SIR, subYears,\n",
    "                        rows_cols_env[0], rows_cols_env[1], rows_cols_env[2], rows_cols_env[3])\n",
    "\n",
    "# Information passed back from \"read_Tb_whole\" reader includes:\n",
    "CETB_SIR = data_SIR['TB']   # 3-D Tb time-series array of TB\n",
    "cal_date = data_SIR['cal_date'] # 1-D array of dates, these will get passed to later functions\n",
    "cal_year = data_SIR['cal_year']    # 1-D array of years\n",
    "cal_month = data_SIR['cal_month']   # 1-D array of months\n",
    "# data_SIR['latitude'], data_SIR['longitude'] # 2-D arrays of subset pixel lat/lons\n",
    "# data_SIR['x'], data_SIR['y'] # 2-D arrays of subset pixel projected x/y\n",
    "# data_SIR['gpd'] # name of EASE2 projection that the subset was derived from\n",
    "\n",
    "# load GRD Tb data\n",
    "data_GRD = read_Tb_whole(datadir_GRD, prefix_GRD, subYears,\n",
    "                         rows_cols_GRD[0], rows_cols_GRD[1], rows_cols_GRD[2], rows_cols_GRD[3])\n",
    "\n",
    "CETB_GRD = data_GRD['TB']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate DAV for the Tb data that was imported\n",
    "DAV_SIR = calc_DAV(data_SIR['TB'])\n",
    "DAV_GRD = calc_DAV(data_GRD['TB'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate histogram - SIR - all data\n",
    "# This histogram will include all points that are in the CETB_SIR array\n",
    "year=2010\n",
    "Tb_threshold=225  #set above by sensor\n",
    "data = CETB_SIR[cal_year==year]\n",
    "#data = CETB_SIR[:] # CETB_data for all pixels in subset\n",
    "data = data[data>=0]\n",
    "bins = range(150,300)  # bins for histogram\n",
    "fig,ax = plt.subplots()\n",
    "ax.hist(data, bins)\n",
    "ax.set_title('TB SIR Histogram, area=%s, year=%s' % (Site, str(year)))\n",
    "ax.axvline(x=Tb_threshold, color='red')\n",
    "ax.set_xlabel('Brightness Temp (K)') \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## create MOD in data frames with geolocation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get array of average MODs for SIR pixels for all the years loaded\n",
    "#Colorado (Johnson et al 2020) used 5 times in 7 day window, for locations in Colorado\n",
    "#Patagonia (Monahan and Ramage 2010) and Yukon (Semmens et al 2013?) used 3 times in 5 day window\n",
    "\n",
    "# From Matias Fall 2022 forward facing indexer to assign rolling sum value to the beginning of window\n",
    "# Setting the window_size sets the number of observations,'14' would be 7 days (2 measurements per day)\n",
    "# If you don't want to use the forward facing indexer, then change \"window\" to a numeral, this will assign\n",
    "# the rolling sum value to the end of the window\n",
    "MOD_window = 10\n",
    "MOD_count = 1\n",
    "indexer = pd.api.indexers.FixedForwardWindowIndexer(window_size=MOD_window)\n",
    "window = indexer   \n",
    "\n",
    "# number of Tb/DAV exceedances to trigger MOD\n",
    "#From Johnson et al 2020 AMSRE rSIR Tb >= 249 DAV>=13 and AMSRE GRD Tb>=243 DAV>=14\n",
    "#From Johnson et al 2020 SSMI rSIR and GRD Tb>=247 DAV>=10\n",
    "DAV_threshold = 10\n",
    "#Tb_threshold = 225 #comment back in for localized testing\n",
    "\n",
    "# Number of Tb/DAV exceedances to trigger EHD = end of high DAV\n",
    "# At the current time, EHD is not using the forward facing indexer, using default behavior\n",
    "# From Matias Fall 2022\n",
    "EHD_window = 20 \n",
    "EHD_count = 7 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# sir MOD array - MOD will be in day of year (DOY) #changeback to Years for all years\n",
    "MOD_DOY_df, meltflag_df, EHD_DOY_df, EHDflag_df = MOD_array(\n",
    "    datadir_SIR, prefix_SIR, data_SIR, DAV_SIR, rows_cols_env, \n",
    "    subYears, window, MOD_count, EHD_window, EHD_count, DAV_threshold, Tb_threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## MOD of the GRD pixel - avg all years\n",
    "MOD_DOY_GRD_df, meltflag_GRD_df, EHD_DOY_GRD_df, EHDflag_GRD_df = MOD_array(\n",
    "    datadir_GRD, prefix_GRD, data_GRD, DAV_GRD, rows_cols_GRD, \n",
    "    subYears, window, MOD_count, EHD_window, EHD_count, DAV_threshold, Tb_threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MOD_DOY_GRD_df "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "EHD_DOY_GRD_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "meltflag_GRD_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EHDflag_GRD_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notes about changes from original notebooks:\n",
    "\n",
    "The old notebooks used to call MOD_array to get the average MOD for a set of years \n",
    "and then call MOD_array_year for a given year of interest.\n",
    "\n",
    "Now we just call MOD_array for SIR data and GRD data, and get back a data frame with \n",
    "MOD columns for each individual year, and one column for the avg MOD for all the years.\n",
    "\n",
    "This will run much faster, and can be saved and just re-read from a saved file on disk.\n",
    "\n",
    "N.B. With this approach, we can create pickle files once, and then \n",
    "read saved pickle files with this command:\n",
    "\n",
    "new = pd.read_pickle(MOD_DOY_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the MOD by year data frames for SIR and GRD to pickle files (with lat/lon geolocation)\n",
    "\n",
    "Also saving geolocation and melt onset flag data frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modsDir = \"%s/MODs\" % outDir\n",
    "if not os.path.isdir(modsDir):\n",
    "    os.makedirs(modsDir)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.MOD.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_SIR['gpd'], region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, \n",
    "    subYears[0], subYears[-1], MOD_count, MOD_window)\n",
    "MOD_DOY_df.to_pickle(filename)\n",
    "print(\"MOD_DOY dataframe saved to %s\\n\" % filename)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.MOD.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_GRD['gpd'], region, SiteLabel, sat_GRD, channel_GRD, alg_GRD, \n",
    "    subYears[0], subYears[-1], MOD_count, MOD_window)\n",
    "MOD_DOY_GRD_df.to_pickle(filename)\n",
    "print(\"MOD_DOY_GRD dataframe saved to %s\\n\" % filename)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.meltflag.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_SIR['gpd'], region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, \n",
    "    subYears[0], subYears[-1], MOD_count, MOD_window)\n",
    "meltflag_df.to_pickle(filename)\n",
    "print(\"meltflag_df dataframe saved to %s\\n\" % filename)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.meltflag.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_GRD['gpd'], region, SiteLabel, sat_GRD, channel_GRD, alg_GRD, \n",
    "    subYears[0], subYears[-1], MOD_count, MOD_window)\n",
    "meltflag_GRD_df.to_pickle(filename)\n",
    "print(\"meltflag_GRD_df dataframe saved to %s\\n\" % filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the MOD by year as geotiff images\n",
    "\n",
    "This will save one geotiff for each year plus one each for the Avg, Min and Max MODs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outbasename = \"%s/%s-%s-%s.%s.%s.%s.MOD.%1dof%02d\" % (\n",
    "    modsDir, data_SIR['gpd'], region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, MOD_count, MOD_window)\n",
    "out = write_df_to_geotiff(MOD_DOY_df, data_SIR['gpd'], outbasename, verbose=True)\n",
    "\n",
    "outbasename = \"%s/%s-%s-%s.%s.%s.%s.MOD.%1dof%02d\" % (\n",
    "    modsDir, data_GRD['gpd'], region, SiteLabel, sat_GRD, channel_GRD, alg_GRD, MOD_count, MOD_window)\n",
    "out = write_df_to_geotiff(MOD_DOY_GRD_df, data_GRD['gpd'], outbasename, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the EHD by year data frames for SIR and GRD to pickle files (with lat/lon geolocation)\n",
    "\n",
    "Also saving geolocation and end of high DAV flag data frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modsDir = \"%s/MODs\" % outDir\n",
    "if not os.path.isdir(modsDir):\n",
    "    os.makedirs(modsDir)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.EHD.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_SIR['gpd'], region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, \n",
    "    subYears[0], subYears[-1], EHD_count, EHD_window)\n",
    "\n",
    "EHD_DOY_df.to_pickle(filename)\n",
    "print(\"EHD_DOY dataframe saved to %s\\n\" % filename)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.EHD.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_GRD['gpd'], region, SiteLabel, sat_GRD, channel_GRD, alg_GRD, \n",
    "    subYears[0], subYears[-1], EHD_count, EHD_window)\n",
    "EHD_DOY_GRD_df.to_pickle(filename)\n",
    "print(\"EHD_DOY_GRD dataframe saved to %s\\n\" % filename)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.EHDflag.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_SIR['gpd'], region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, \n",
    "    subYears[0], subYears[-1], EHD_count, EHD_window)\n",
    "EHDflag_df.to_pickle(filename)\n",
    "print(\"meltflag_df dataframe saved to %s\\n\" % filename)\n",
    "\n",
    "filename = \"%s/%s-%s-%s.%s.%s.%s.%s-%s.EHDflag.%1dof%02d.pkl\" % (\n",
    "    modsDir, data_GRD['gpd'], region, SiteLabel, sat_GRD, channel_GRD, alg_GRD, \n",
    "    subYears[0], subYears[-1], EHD_count, EHD_window)\n",
    "EHDflag_GRD_df.to_pickle(filename)\n",
    "print(\"EHDflag_GRD_df dataframe saved to %s\\n\" % filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the EHD by year as geotiff images\n",
    "\n",
    "This will save one geotiff for each year plus one each for the Avg, Min and Max EHDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outbasename = \"%s/%s-%s-%s.%s.%s.%s.EHD.%1dof%02d\" % (\n",
    "    modsDir, data_SIR['gpd'], region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, EHD_count, EHD_window)\n",
    "out = write_df_to_geotiff(EHD_DOY_df, data_SIR['gpd'], outbasename, verbose=True)\n",
    "\n",
    "outbasename = \"%s/%s-%s-%s.%s.%s.%s.EHD.%1dof%02d\" % (\n",
    "    modsDir, data_GRD['gpd'], region, SiteLabel, sat_GRD, channel_GRD, alg_GRD, EHD_count, EHD_window)\n",
    "out = write_df_to_geotiff(EHD_DOY_GRD_df, data_GRD['gpd'], outbasename, verbose=True)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Set your parameters for the pkl file to read\n",
    "\n",
    "subYears = Years[0:2]\n",
    "#subYears = Years\n",
    "\n",
    "# a short string in the filename to indicate whether window was the indexer function or a plain integer\n",
    "#MODrollingSumLabel = 'end'\n",
    "MODrollingSumLabel = 'beg'\n",
    "EHDrollingSumLabel = 'end'\n",
    "#EHDrollingSumLabel = 'beg'\n",
    "\n",
    "MOD_window = 10\n",
    "MOD_count = 3\n",
    "DAV_threshold = 10\n",
    "Tb_threshold = 225\n",
    "EHD_window = 20 \n",
    "EHD_count = 7 \n",
    "    \n",
    "MODinfo = \"MOD_C%1dW%02d%sT%03dD%02d\" % (\n",
    "    MOD_count, MOD_window, MODrollingSumLabel, Tb_threshold, DAV_threshold) \n",
    "EHDinfo = \"EHD_C%1dW%02d%sT%03dD%02d\" % (\n",
    "    EHD_count, EHD_window, EHDrollingSumLabel, Tb_threshold, DAV_threshold)\n",
    "\n",
    "meltflaginfo = \"meltflag_C%1dW%02d%sT%03dD%02d\" % (\n",
    "    MOD_count, MOD_window, MODrollingSumLabel, Tb_threshold, DAV_threshold) \n",
    "EHDflaginfo = \"EHDflag_C%1dW%02d%sT%03dD%02d\" % (\n",
    "    EHD_count, EHD_window, EHDrollingSumLabel, Tb_threshold, DAV_threshold) \n",
    "\n",
    "sirMODfile = \"%s/%s.%s.%s.%s.%s.%s-%s.%s.pkl\" % (\n",
    "    pklDir, gpd_SIR, region, SiteLabel, sat_SIR, channel_SIR, \n",
    "    subYears[0], subYears[-1], MODinfo)\n",
    "grdMODfile = \"%s/%s.%s.%s.%s.%s.%s-%s.%s.pkl\" % (\n",
    "    pklDir, gpd_GRD, region, SiteLabel, sat_GRD, channel_GRD, \n",
    "    subYears[0], subYears[-1], MODinfo)\n",
    "\n",
    "sirmeltflagfile = \"%s/%s.%s.%s.%s.%s.%s-%s.%s.pkl\" % (\n",
    "    pklDir, gpd_SIR, region, SiteLabel, sat_SIR, channel_SIR, \n",
    "    subYears[0], subYears[-1], meltflaginfo)\n",
    "grdmeltflagfile = \"%s/%s.%s.%s.%s.%s.%s-%s.%s.pkl\" % (\n",
    "    pklDir, gpd_GRD, region, SiteLabel, sat_GRD, channel_GRD, \n",
    "    subYears[0], subYears[-1], meltflaginfo)\n",
    "\n",
    "MOD_DOY_GRD_df = pd.read_pickle(grdMODfile)\n",
    "meltflag_GRD_df = pd.read_pickle(grdmeltflagfile)\n",
    "\n",
    "MOD_DOY_df = pd.read_pickle(sirMODfile)\n",
    "meltflag_df = pd.read_pickle(sirmeltflagfile)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display maps of the MOD data\n",
    "\n",
    "This next example displays the requested SIR MOD map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#how to create subplots with basemap (w maps) https://basemaptutorial.readthedocs.io/en/latest/subplots.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLOT the MOD MAP - avg of all years or one year of interest - SIR\n",
    "# Choose plot type - 'year' for one year of interest, and specify the year,\n",
    "# or 'all' for the average of all years\n",
    "plot_type='year'  \n",
    "year=2010\n",
    "#plot_type='all'\n",
    "\n",
    "if plot_type == 'all':\n",
    "    array = MOD_DOY_df['Avg']\n",
    "    label = 'AvgDOY'+str(subYears[0])+'-'+str(subYears[-1])\n",
    "    title = sensor_SIR + '-' + sat_SIR + '-' + Site + ' - MOD - (Avg DOY '+str(subYears[0])+'-'+str(subYears[-1])+')'\n",
    "elif plot_type == 'year':\n",
    "    array = MOD_DOY_df[year]\n",
    "    label = str(year)+'DOY' \n",
    "    title = sensor_SIR + '-' +sat_SIR + '-'+ Site + ' - MOD - (' + str(year) + ' DOY)'\n",
    "\n",
    "# Set a few common things\n",
    "graticule_fontsize = 4\n",
    "\n",
    "# create figure and axes instances\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "ax = fig.add_axes([0.1,0.1,0.8,0.8])\n",
    "#ax = fig.add_subplot(121)\n",
    "# create polar stereographic Basemap instance.\n",
    "m = Basemap(projection='stere',\n",
    "            lon_0=(lon_start+lon_end)/2,\n",
    "            lat_0=(lat_start+lat_end)/2,\n",
    "            lat_ts=90.,\n",
    "            llcrnrlat=(lat_start-.3),\n",
    "            urcrnrlat=(lat_end+.5),\n",
    "            llcrnrlon=(lon_start-.3),\n",
    "            urcrnrlon=(lon_end+.3),\n",
    "            rsphere=6371200.,\n",
    "            resolution='l',\n",
    "            area_thresh=10000, \n",
    "            epsg=3857)\n",
    "m.drawcoastlines()\n",
    "m.drawcountries()\n",
    "m.drawstates()\n",
    "#m.etopo(scale=5, alpha=0.7)\n",
    "#m.bluemarble(scale=1)\n",
    "parallels = np.arange(0.,90,0.5)\n",
    "m.drawparallels(parallels, labels=[1,0,0,0], fontsize=graticule_fontsize)\n",
    "# draw meridians\n",
    "meridians = np.arange(0,360.,0.5)\n",
    "m.drawmeridians(meridians, labels=[0,0,0,1], fontsize=graticule_fontsize)\n",
    "m.arcgisimage(service='World_Physical_Map', xpixels = 1500, verbose= True)\n",
    "\n",
    "#specific shapefile overlays can be added here:\n",
    "#m.readshapefile('/home/mij216/ExploringCETB/IN_Hunza_at_DainyorBridge', 'Hunza', color='red', linewidth=1)\n",
    "#m.readshapefile('/mnt/data3/rgi6.0/06_rgi60_Iceland', 'Iceland', color='red', linewidth=1)\n",
    "\n",
    "x,y = m(MOD_DOY_df['longitude'].values, MOD_DOY_df['latitude'].values)\n",
    "# scatter plot size argument is size in points**2, default\n",
    "# is rcParams['lines.markersize'] ** 2.\n",
    "scatter_size = 24.\n",
    "im = m.scatter(x, y, c=array, s=scatter_size, marker='s',\n",
    "              lw=0, cmap='BuPu_r', alpha=.6) \n",
    "plt.title(title)\n",
    "\n",
    "# Get the current axes and configure placement so colorbar will\n",
    "# be 5% of ax width and the padding will be 0.05 in.\n",
    "ax = plt.gca()\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax, label='DOY')\n",
    "plt.clim(75,175)  #color bar limits\n",
    "\n",
    "modsDir = \"%s/MODs\" % outDir\n",
    "if not os.path.isdir(modsDir):\n",
    "    os.makedirs(modsDir)\n",
    "filename = \"%s/%s.%s.%s.%s.%s.%s.png\" % (\n",
    "    modsDir, region, SiteLabel, sat_SIR, channel_SIR, alg_SIR, label)\n",
    "plt.savefig(filename)\n",
    "print(\"image saved to %s\\n\" % filename)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#color bars here https://matplotlib.org/stable/tutorials/colors/colormaps.html\n",
    "#documentation on basemap https://matplotlib.org/basemap/users/laea.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map of GRD (left) and SIR (right) MOD maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Trial-and-error decide on size for GRD vs. SIR scatter markers\n",
    "# scatter size is in points**2, with default of\n",
    "# rcParams['lines.markersize'] ** 2.\n",
    "#np.sqrt((1200)/8)^2\n",
    "sGRD=700\n",
    "sSIR=12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2plots left (GRD) and right (SIR)\n",
    "#redo basemap to be in easegrid projection EASE2.0\n",
    "#\n",
    "# Choose plot type - 'year' for one year of interest, and specify the year,\n",
    "# or 'all' for the average of all years\n",
    "plot_type='year'  \n",
    "year=2010\n",
    "#plot_type='all'\n",
    "\n",
    "\n",
    "if plot_type == 'all':\n",
    "    label = 'MOD.Avg.%1dof%02d.%4d-%4d' % (MOD_count, MOD_window, subYears[0], subYears[-1])\n",
    "    sir_array = MOD_DOY_df['Avg']\n",
    "    sir_title = sensor_SIR + '-' + sat_SIR + '-' + Site + ' rSIR MOD (Avg DOY '+str(subYears[0])+'-'+str(subYears[-1])+')'\n",
    "    grd_array = MOD_DOY_GRD_df['Avg']\n",
    "    grd_title = sensor_GRD + '-' + sat_GRD + '-' + Site + ' GRD MOD (Avg DOY '+str(subYears[0])+'-'+str(subYears[-1])+')'\n",
    "elif plot_type == 'year':\n",
    "    label = 'MOD.%1dof%02d.%4d' % (MOD_count, MOD_window, year)\n",
    "    sir_array = MOD_DOY_df[year]\n",
    "    sir_title = sensor_SIR + '-' + sat_SIR + '-' + Site + ' rSIR MOD (' + str(year) + ' DOY)'\n",
    "    grd_array = MOD_DOY_GRD_df[year]\n",
    "    grd_title = sensor_GRD + '-' + sat_GRD + '-' + Site + ' GRD MOD (' + str(year) + ' DOY)'\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14,5))\n",
    "\n",
    "# Set a few common things\n",
    "graticule_fontsize = 7\n",
    "\n",
    "#min day of year to plot\n",
    "minday=75\n",
    "#max day of year to plot\n",
    "maxday=175\n",
    "\n",
    "# create polar stereographic Basemap instance.\n",
    "axes[0].set_title(grd_title)\n",
    "m = Basemap(ax=axes[0],projection='stere',\n",
    "            lon_0=(lon_start+lon_end)/2,\n",
    "            lat_0=(lat_start+lat_end)/2,lat_ts=90.,\n",
    "            llcrnrlat=(lat_start-.3),\n",
    "            urcrnrlat=(lat_end+.5),\n",
    "            llcrnrlon=(lon_start-.3),\n",
    "            urcrnrlon=(lon_end+.3),\n",
    "            rsphere=6371200.,\n",
    "            resolution='l',\n",
    "            area_thresh=10000, \n",
    "            epsg=3857)\n",
    "m.drawcoastlines()\n",
    "m.drawcountries()\n",
    "#m.drawstates()\n",
    "#m.etopo(scale=3, alpha=0.7)\n",
    "#m.bluemarble(scale=1)\n",
    "parallels = np.arange(0.,90,1.0)\n",
    "m.drawparallels(parallels, labels=[1,0,0,0], fontsize=graticule_fontsize)\n",
    "# draw meridians\n",
    "meridians = np.arange(0,360.,1.0)\n",
    "m.drawmeridians(meridians, labels=[0,0,0,1], fontsize=graticule_fontsize)\n",
    "m.arcgisimage(service='World_Physical_Map', xpixels = 1500, verbose= True)\n",
    "\n",
    "#m.readshapefile('/home/mij216/ExploringCETB/IN_Hunza_at_DainyorBridge', 'Hunza', color='red', linewidth=1)\n",
    "#m.readshapefile('/mnt/data3/rgi6.0/06_rgi60_Iceland', 'Vatnajokull', color='red', linewidth=1)\n",
    "\n",
    "x,y = m(MOD_DOY_GRD_df['longitude'].values, MOD_DOY_GRD_df['latitude'].values)\n",
    "im0 = m.scatter(x, y, c=grd_array, \n",
    "                s=sGRD, marker='s',lw=0,cmap='BuPu_r', alpha=.6, vmin=minday, vmax=maxday) \n",
    "\n",
    "# Get the current axes and configure placement so colorbar will\n",
    "# be 5% of ax width and the padding will be 0.08 in.\n",
    "divider = make_axes_locatable(axes[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.08)\n",
    "cbar = plt.colorbar(im0, cax=cax, label='DOY')\n",
    "\n",
    "axes[1].set_title(sir_title)\n",
    "m = Basemap(ax=axes[1],projection='stere',\n",
    "            lon_0=(lon_start+lon_end)/2,\n",
    "            lat_0=(lat_start+lat_end)/2,lat_ts=90.,\n",
    "            llcrnrlat=(lat_start-.3),\n",
    "            urcrnrlat=(lat_end+.5),\n",
    "            llcrnrlon=(lon_start-.3),\n",
    "            urcrnrlon=(lon_end+.3),\n",
    "            rsphere=6371200.,\n",
    "            resolution='l',\n",
    "            area_thresh=10000, \n",
    "            epsg=3857)\n",
    "m.drawcoastlines()\n",
    "m.drawcountries()\n",
    "#m.drawstates()\n",
    "#m.etopo(scale=3, alpha=0.7)\n",
    "#m.bluemarble(scale=1)\n",
    "parallels = np.arange(0.,90,1.0)\n",
    "m.drawparallels(parallels, labels=[1,0,0,0], fontsize=graticule_fontsize)\n",
    "# draw meridians\n",
    "meridians = np.arange(0,360.,1.0)\n",
    "m.drawmeridians(meridians, labels=[0,0,0,1], fontsize=graticule_fontsize)\n",
    "m.arcgisimage(service='World_Physical_Map', xpixels = 1500, verbose= True)\n",
    "\n",
    "#m.readshapefile('/home/mij216/ExploringCETB/IN_Hunza_at_DainyorBridge', 'Hunza', color='red', linewidth=1)\n",
    "#m.readshapefile('/mnt/data3/rgi6.0/06_rgi60_Iceland', 'Vatnajokull', color='red', linewidth=1)\n",
    "\n",
    "x,y = m(MOD_DOY_df['longitude'].values, MOD_DOY_df['latitude'].values)\n",
    "im1 = m.scatter(x, y, c=sir_array, \n",
    "                s=sSIR, marker='s', lw=0,cmap='BuPu_r', alpha=.6,vmin=minday,vmax=maxday) \n",
    "#fig.colorbar(im1, ax=axes[1],label='DOY')\n",
    "#plt.clim(75,175)  #color bar limits\n",
    "\n",
    "# Get the current axes and configure placement so colorbar will\n",
    "# be 5% of ax width and the padding will be 0.05 in.\n",
    "divider = make_axes_locatable(axes[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "cbar = plt.colorbar(im1, cax=cax, label='DOY')\n",
    "\n",
    "plt.tight_layout(pad=3.0)\n",
    "\n",
    "\n",
    "filename = \"%s/%s.%s.%s.%s.%s.%s.png\" % (\n",
    "    modsDir, region, SiteLabel, sat_SIR, channel_SIR, \"GRDvsSIR\", label)\n",
    "plt.savefig(filename)\n",
    "print(\"image saved to %s\\n\" % filename)\n",
    "\n",
    "# Apparently order is important here: you have to save the figure before you do plt.show()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#setup a dataframe so that all the data from ssmi and ssmis are read into one array for long term average\n",
    "#need to decide which sensor has priority when there is overlap\n",
    "#can create a reader that loops through enough sensors to get a full time series \n",
    "#(skip F10- orbit very elliptical and F19 - short)\n",
    "#create that time series for the analysis, save for analysis\n",
    "#\n",
    "#save dataframes with MOD so that we can do other analysys -- done\n",
    "#pickle the data frame that is inside the MOD_array_DOY function -- done \n",
    "#\n",
    "#check in CETB_analysis.py that when it does the count/window it is in fact \n",
    "#saving the first occurrence of the count, even though it doesn't know until it gets to 3x  -- pretty sure this is done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Add to this\n",
    "#Plot of which pixels have melted as of a specific day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure out refreeze? then from MOD to Freeze, count # of dates with Tb>=threshold (with or without DAV? might depend on site)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Add to this\n",
    "#Plot of how many days (or occurrences) experienced melt -- this would need to read the meltflag_df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Add to this\n",
    "#Plot of how many short (1-2 occurrences) melt days EMEs prior to the first MOD\n",
    "#From Jan 1 - MOD, Number of days with Tb>= threshold and DAV>= threshold (same thresholds) --this would come from meltflag_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}