Modeled - Tmax and Tmin - 20 year - monthly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cd9e6c27",
   "metadata": {},
   "source": [
    "## WICCI Downscaled Data: Tmax and Tmin - 20 year - monthly\n",
    "\n",
    "Objectives\n",
    "* aggregate data for tmax and tmin, e.g. mean, standard deviation, days over 100, for RCP4.5 and RCP8.5\n",
    "* run for every 20 year model window, e.g. 2021-2040 for 2030 average\n",
    "* calculate averages by whole month increments\n",
    "* creat new netcdf file(s) for aggregate data for each 20-year timeframe\n",
    "* calculate intermodel standard deviation differently -- find average of each model within the current time window and\n",
    "  then calculate standard deviation across models\n",
    "\n",
    "Eric Compas, compase@uww.edu 11/17/2021, 1/19/2022, 1/25/2022, 2/28/2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b488720d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import netCDF4\n",
    "import numpy as np\n",
    "import os\n",
    "import datetime\n",
    "import gc\n",
    "from netCDF4 import Dataset,num2date,date2num"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3b7aaae1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae1ccaba",
   "metadata": {},
   "source": [
    "Which folder are we getting the downscaled climate data from? Needs to be the data from WICCI Climate Working Group with folder structure as provided. This level should contain the 'rcp45' and 'rcp85' folders to work correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ff913133",
   "metadata": {},
   "outputs": [],
   "source": [
    "#base_folder = \"Z:/Climate_Data\"\n",
    "base_folder = \"F:/Climate_Data\"\n",
    "#base_folder = \"/Users/ericcompas/Climate_Data\"\n",
    "if not os.path.isdir(base_folder):\n",
    "    print(\"Base folder not valid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0384094e",
   "metadata": {},
   "source": [
    "Folder to store the resulting netCDF files in. Can be anywhere on your file system. Note that this script will not overwrite existing files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c3853f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "#out_folder = \"Z:/Climate_Data/aggregate_files\"\n",
    "out_folder = \"C:/Users/Eric/Dropbox/Climate_Data/aggregate_files\"\n",
    "#out_folder  = \"/Users/ericcompas/Climate_Data/output\"\n",
    "if not os.path.isdir(out_folder):\n",
    "    print(\"Out folder not valid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76f0f3e",
   "metadata": {},
   "source": [
    "The BIG LOOP.\n",
    "\n",
    "Look across all climate scenarios, years, models, and create aggregate file for each set of years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f6d6f533",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing climate scenario rcp85\n",
      "  Processing base year 2030\n",
      "    Processing GCM ACCESS1-0\n",
      "      Processing year 2021\n",
      "      Processing year 2022\n",
      "      Processing year 2023\n",
      "      Processing year 2024\n",
      "      Processing year 2025\n",
      "      Processing year 2026\n",
      "      Processing year 2027\n",
      "      Processing year 2028\n",
      "      Processing year 2029\n",
      "      Processing year 2030\n",
      "      Processing year 2031\n",
      "      Processing year 2032\n",
      "      Processing year 2033\n",
      "      Processing year 2034\n",
      "      Processing year 2035\n",
      "      Processing year 2036\n",
      "      Processing year 2037\n",
      "      Processing year 2038\n",
      "      Processing year 2039\n",
      "      Processing year 2040\n",
      "    Processing GCM ACCESS1-3\n",
      "      Processing year 2021\n",
      "      Processing year 2022\n",
      "      Processing year 2023\n",
      "      Processing year 2024\n",
      "      Processing year 2025\n",
      "      Processing year 2026\n",
      "      Processing year 2027\n",
      "      Processing year 2028\n",
      "      Processing year 2029\n",
      "      Processing year 2030\n",
      "      Processing year 2031\n",
      "      Processing year 2032\n",
      "      Processing year 2033\n",
      "      Processing year 2034\n",
      "      Processing year 2035\n",
      "      Processing year 2036\n",
      "      Processing year 2037\n",
      "      Processing year 2038\n",
      "      Processing year 2039\n",
      "      Processing year 2040\n",
      "    Processing GCM CanESM2\n",
      "      Processing year 2021\n",
      "      Processing year 2022\n",
      "      Processing year 2023\n",
      "      Processing year 2024\n",
      "      Processing year 2025\n",
      "      Processing year 2026\n",
      "      Processing year 2027\n",
      "      Processing year 2028\n",
      "      Processing year 2029\n",
      "      Processing year 2030\n",
      "      Processing year 2031\n",
      "      Processing year 2032\n",
      "      Processing year 2033\n",
      "      Processing year 2034\n",
      "      Processing year 2035\n",
      "      Processing year 2036\n",
      "      Processing year 2037\n",
      "      Processing year 2038\n",
      "      Processing year 2039\n",
      "      Processing year 2040\n",
      "Calc mean for index start 0 to 31\n",
      "Calc mean for index start 31 to 59\n",
      "Calc mean for index start 59 to 90\n",
      "Calc mean for index start 90 to 120\n",
      "Calc mean for index start 120 to 151\n",
      "Calc mean for index start 151 to 181\n",
      "Calc mean for index start 181 to 212\n",
      "Calc mean for index start 212 to 243\n",
      "Calc mean for index start 243 to 273\n",
      "Calc mean for index start 273 to 304\n",
      "Calc mean for index start 304 to 334\n",
      "Calc mean for index start 334 to 366\n"
     ]
    }
   ],
   "source": [
    "# loop across models\n",
    "#models = [\"rcp45\",\"rcp85\"]\n",
    "models = [\"rcp85\"]\n",
    "for m in models:\n",
    "    print(\"Processing climate scenario \"+m)\n",
    "\n",
    "    # get subfolders/GCMs for current model\n",
    "    GCMs = os.listdir(os.path.join(base_folder,m))\n",
    "    GCMs = GCMs[:3]  # for testing only\n",
    "\n",
    "    # loop across years in 20-year increments\n",
    "    for y in range(2030,2091,20):    # 2030,2091,20\n",
    "        print(\"  Processing base year \"+str(y))\n",
    "\n",
    "        # create (large) arrays for all data within current set of years\n",
    "        # loop through GCMS first and then years in current time window\n",
    "        # (create large blank masked arrays more memory efficient?)\n",
    "        num_models = len(GCMs)\n",
    "        years = 20\n",
    "        realizations = 3\n",
    "        tmax_all_array = np.ma.empty([num_models*years*realizations, 365, 48, 62])\n",
    "        tmax_means_array = np.ma.empty([num_models, 365, 48, 62])\n",
    "        tmin_all_array = np.ma.empty([num_models*years*realizations, 365, 48, 62])\n",
    "        tmin_means_array = np.ma.empty([num_models, 365, 48, 62])\n",
    "        \n",
    "        # indices for adding to blank arrays\n",
    "        i = 0\n",
    "        j = 0\n",
    "        \n",
    "        for gcm in GCMs:    \n",
    "            print(\"    Processing GCM \"+gcm)\n",
    "            tmax_gcm_array = []\n",
    "            tmin_gcm_array = []\n",
    "            for yr in range(y-9,y+11):\n",
    "                print(\"      Processing year \"+str(yr))\n",
    "                realizations = [\"01\",\"02\",\"03\"]\n",
    "                # get data for all three \"realizations\" from each GCM\n",
    "                for r in realizations:\n",
    "                    nf = os.path.join(base_folder,m,gcm,\"r1i1p1\",\"temp_\"+r+\"_\"+str(yr)+\".nc\")\n",
    "                    try:\n",
    "                        n = netCDF4.Dataset(nf)\n",
    "                    except:\n",
    "                        # only happens if file is missing -- does happen for a couple of years in the \n",
    "                        # first 20-year window\n",
    "                        q = 0    # does nothing\n",
    "                    else:\n",
    "                        # found file, now lets get tmin and tmax for this GCM realization\n",
    "                        tmax = n.variables['tmax']\n",
    "                        tmin = n.variables['tmin']\n",
    "                        \n",
    "                        # convert to F\n",
    "                        tmax = (tmax[:,:,:] * 1.8) + 32\n",
    "                        tmin = (tmin[:,:,:] * 1.8) + 32\n",
    "                        \n",
    "                        # get mask - to reapply below\n",
    "                        mask = np.ma.getmask(tmax[0])\n",
    "                        \n",
    "                        # handle leap years -- average Feb 28 and 29\n",
    "                        if tmax.shape[0] == 366:\n",
    "                            # get and average Feb 28 and 29 for tmax\n",
    "                            tmax_leap = tmax[58:60,:,:]\n",
    "                            tmax_mean_leap = np.ma.mean(tmax_leap,axis=0)\n",
    "                            tmax_part1 = np.ma.append(tmax[0:58,:,:],np.ma.expand_dims(tmax_mean_leap,axis=0),axis=0)\n",
    "                            tmax = np.ma.append(tmax_part1,tmax[60:367,:,:],axis=0)\n",
    "                            # ...and for tmin\n",
    "                            tmin_leap = tmin[58:60,:,:]\n",
    "                            tmin_mean_leap = np.ma.mean(tmin_leap,axis=0)\n",
    "                            tmin_part1 = np.ma.append(tmin[0:58,:,:],np.ma.expand_dims(tmin_mean_leap,axis=0),axis=0)\n",
    "                            tmin = np.ma.append(tmin_part1,tmin[60:367,:,:],axis=0)\n",
    "                        \n",
    "                        # add tmax array for gcm, year, realization to total array\n",
    "                        tmax_all_array[i] = tmax[:,:,:]\n",
    "                        tmin_all_array[i]= tmin[:,:,:]\n",
    "                        \n",
    "                        # add tmin and tmax array for gcm, year, realization to gcm-specific array\n",
    "                        # (used below for inter-model standard deviation)\n",
    "                        tmax_gcm_array.append(tmax[:,:,:])\n",
    "                        tmin_gcm_array.append(tmin[:,:,:])\n",
    "                        \n",
    "                        i += 1\n",
    "                        \n",
    "                    # end realization loop\n",
    "                # end year loop\n",
    "            \n",
    "            # calculate mean temperature for gcm over time window (for std calculation)\n",
    "            tmax_model_mean = np.ma.mean(tmax_gcm_array,axis=0)\n",
    "            tmin_model_mean = np.ma.mean(tmin_gcm_array,axis=0)\n",
    "\n",
    "            # add this to an array of means (one for each gcm)\n",
    "            tmax_means_array[j] = tmax_model_mean[:,:,:]\n",
    "            tmin_means_array[j] = tmin_model_mean[:,:,:]\n",
    "            \n",
    "            j += 1\n",
    "            \n",
    "            # end GCM loop\n",
    "            \n",
    "        ##################################################\n",
    "        ## calculate metrics for current 20-year window ##\n",
    "        ##################################################\n",
    "        \n",
    "        # calculate daily mean from total array    \n",
    "        tmax_mean = np.ma.mean(tmax_all_array,axis=0)\n",
    "        tmin_mean = np.ma.mean(tmin_all_array,axis=0)\n",
    "        \n",
    "        # convert/calculate monthly values from total array (in days)\n",
    "        tmax_months_mean = np.ma.empty([12, 48, 62])\n",
    "        tmin_months_mean = np.ma.empty([12, 48, 62])\n",
    "        tmax_months_std = np.ma.empty([12, 48, 62])\n",
    "        tmin_months_std = np.ma.empty([12, 48, 62])\n",
    "        tmax_months_gt100 = np.ma.empty([12, 48, 62])\n",
    "        tmax_months_gt95 = np.ma.empty([12, 48, 62])\n",
    "        tmax_months_gt90 = np.ma.empty([12, 48, 62])\n",
    "        tmax_months_comfort = np.ma.empty([12, 48, 62])\n",
    "        tmin_months_gt70 = np.ma.empty([12, 48, 62])\n",
    "        tmin_months_lt32 = np.ma.empty([12, 48, 62])\n",
    "        tmin_months_lt0 = np.ma.empty([12, 48, 62])\n",
    "        \n",
    "        # index for julian/ordinal calendar (start/end of each month)\n",
    "        julian = [0,31,59,90,120,151,181,212,243,273,304,334,366]\n",
    "        \n",
    "        # loop through months\n",
    "        for month in range(1,13):\n",
    "            # calculate mean values for each month, add to monthly value array\n",
    "            print('Calc mean for index start '+str(julian[month-1])+' to '+str(julian[month]))\n",
    "            tmax_month_mean = np.ma.mean(tmax_mean[(julian[month-1]):(julian[month])],axis=0)\n",
    "            tmax_months_mean[month-1] = tmax_month_mean\n",
    "            tmin_month_mean = np.ma.mean(tmin_mean[(julian[month-1]):(julian[month])],axis=0)\n",
    "            tmin_months_mean[month-1] = tmin_month_mean\n",
    "            \n",
    "            # calculate inter-model standard deviations\n",
    "            all_std_tmax_nparray = np.array(tmax_means_array)\n",
    "            month_std_tmax_nparray = all_std_tmax_nparray[:,(julian[month-1]):(julian[month]),:,:]\n",
    "            all_std_tmin_nparray = np.array(tmin_means_array)\n",
    "            month_std_tmin_nparray = all_std_tmin_nparray[:,(julian[month-1]):(julian[month]),:,:]\n",
    "            # loop over each model -- add to total monthly sample array\n",
    "            for i in range(0, (np.shape(all_std_tmax_nparray)[0])):\n",
    "                model_tmax_array = month_std_tmax_nparray[i]\n",
    "                model_tmin_array = month_std_tmin_nparray[i]\n",
    "                if i==0:\n",
    "                    concat_tmax_array = model_tmax_array\n",
    "                    concat_tmin_array = model_tmin_array\n",
    "                else:\n",
    "                    concat_tmax_array = np.ma.concatenate((concat_tmax_array,model_tmax_array),axis=0)\n",
    "                    concat_tmin_array = np.ma.concatenate((concat_tmin_array,model_tmin_array),axis=0)\n",
    "            tmax_std_month = np.ma.std(concat_tmax_array,axis=0)\n",
    "            tmax_months_std[month-1] = np.ma.masked_array(tmax_std_month,mask)\n",
    "            tmin_std_month = np.ma.std(concat_tmin_array,axis=0)\n",
    "            tmin_months_std[month-1] = np.ma.masked_array(tmin_std_month,mask)\n",
    "            \n",
    "            # calc extreme estimates - tmax\n",
    "            all_nparray = np.array(tmax_all_array)\n",
    "            month_nparray = all_nparray[:,(julian[month-1]):(julian[month]),:,:]\n",
    "            # get number of days in month and number of \"samples\"(needed for estimate formula)\n",
    "            num_days = np.shape(month_nparray)[1]\n",
    "            num_samples = np.shape(month_nparray)[0]*np.shape(month_nparray)[1]     \n",
    "            \n",
    "            # tmax extreme - tmax gt 100 F\n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = (month_nparray >= 100).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmax_months_gt100[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "            \n",
    "            # tmax extreme - tmax gt 95 F     \n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = (month_nparray >= 95).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmax_months_gt95[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "            \n",
    "            # tmax extreme - tmax gt 90 F     \n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = (month_nparray >= 90).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmax_months_gt90[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "            \n",
    "            # tmax extreme(?) - tmax gt 70 F and lt 85 F    \n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = ((month_nparray >= 70) & (month_nparray <= 85)).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmax_months_comfort[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "            \n",
    "            # calc extreme estimates - tmin\n",
    "            all_nparray = np.array(tmin_all_array)\n",
    "            month_nparray = all_nparray[:,(julian[month-1]):(julian[month]),:,:]\n",
    "            # get number of days in month and number of \"samples\"(needed for estimate formula)\n",
    "            num_days = np.shape(month_nparray)[1]\n",
    "            num_samples = np.shape(month_nparray)[0]*np.shape(month_nparray)[1] \n",
    "            \n",
    "            # tmin extreme - tmin gt 70 F     \n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = (month_nparray >= 70).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmin_months_gt70[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "            \n",
    "            # tmin extreme - tmin lt 32 F     \n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = (month_nparray <= 32).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmin_months_lt32[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "            \n",
    "            # tmin extreme - tmin lt 0 F     \n",
    "            # get number of days meeting threshold (by iteration)\n",
    "            extreme_sum_days = (month_nparray <= 0).sum(axis=0)\n",
    "            # add number of days (by month)\n",
    "            extreme_sum = extreme_sum_days.sum(axis=0)\n",
    "            # estimate of number of monthly extreme days\n",
    "            extreme_days = (extreme_sum / num_samples) * num_days\n",
    "            # add to array\n",
    "            tmin_months_lt0[month-1] = np.ma.masked_array(extreme_days,mask)\n",
    "        \n",
    "        break # for testing\n",
    "        \n",
    "        # write netcdf files with results\n",
    "        filename = \"temp_\"+m+\"_\"+str(y)+\"_20yr_monthly.nc\"\n",
    "        newfile = os.path.join(out_folder,filename)\n",
    "        ncfile = netCDF4.Dataset(newfile,mode='w',format='NETCDF4_CLASSIC')\n",
    "        lat_dim = ncfile.createDimension('lat', 48)     # latitude axis\n",
    "        lon_dim = ncfile.createDimension('lon', 62)    # longitude axis\n",
    "        time_dim = ncfile.createDimension('time', None) # unlimited axis (can be appended to).\n",
    "\n",
    "        ncfile.title='Aggregate monthly tmax and tmin values for WICCI downscaled climate data for all GCMs for '+m+' and 20-year window around year '+str(y)\n",
    "        ncfile.subtitle=\"Data source: UW-Madison WICCI; Data aggregation: Eric Compas, compase@uww.edu\"\n",
    "        lat = ncfile.createVariable('lat', np.float64, ('lat',))\n",
    "        lat.units = 'degrees_north'\n",
    "        lat.long_name = 'latitude'\n",
    "        lon = ncfile.createVariable('lon', np.float64, ('lon',))\n",
    "        lon.units = 'degrees_east'\n",
    "        lon.long_name = 'longitude'\n",
    "        time = ncfile.createVariable('time', np.float64, ('time',))\n",
    "        timeunits = 'days since '+str(y)+'-01-01'\n",
    "        time.units = timeunits\n",
    "        time.long_name = 'time'\n",
    "\n",
    "        temp_tmax_mean = ncfile.createVariable('tmax_mean',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmax_mean.units = 'F' # Fahrenheit \n",
    "        temp_tmax_mean.standard_name = 'mean of maximum daily temperature (F) across 20-year window'\n",
    "        temp_tmax_mean.missing_value = -32768\n",
    "        \n",
    "        temp_tmax_std = ncfile.createVariable('tmax_std',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmax_std.units = 'F' # Fahrenheit \n",
    "        temp_tmax_std.standard_name = 'standard deviation of maximum daily temperature (F) across 20-year window'\n",
    "        temp_tmax_std.missing_value = -32768\n",
    "        \n",
    "        temp_tmin_mean = ncfile.createVariable('tmin_mean',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmin_mean.units = 'F' # Fahrenheit \n",
    "        temp_tmin_mean.standard_name = 'mean of minimum daily temperature (F) across 20-year window'\n",
    "        temp_tmin_mean.missing_value = -32768\n",
    "\n",
    "        temp_tmin_std = ncfile.createVariable('tmin_std',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmin_std.units = 'F' # Fahrenheit \n",
    "        temp_tmin_std.standard_name = 'standard deviation of minimum daily temperature (F) across 20-year window'\n",
    "        temp_tmin_std.missing_value = -32768\n",
    "        \n",
    "        temp_tmax_gt100 = ncfile.createVariable('tmax_gt100',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmax_gt100.units = 'days' # number of days\n",
    "        temp_tmax_gt100.standard_name = 'estimated number of days where maximum daily temperature equals or exceeds 100 F '+\\\n",
    "                                        'across 20-year window'\n",
    "        temp_tmax_gt100.missing_value = -32768\n",
    "        \n",
    "        temp_tmax_gt95 = ncfile.createVariable('tmax_gt95',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmax_gt95.units = 'days' # number of days\n",
    "        temp_tmax_gt95.standard_name = 'estimated number of days where maximum daily temperature equals or exceeds 95 F '+ \\\n",
    "                                       'across 20-year window'\n",
    "        temp_tmax_gt95.missing_value = -32768\n",
    "        \n",
    "        temp_tmax_gt90 = ncfile.createVariable('tmax_gt90',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmax_gt90.units = 'days' # number of days\n",
    "        temp_tmax_gt90.standard_name = 'estimated number of days where maximum daily temperature (C) equals or exceeds 90 F '+ \\\n",
    "                                       'across 20-year window'\n",
    "        temp_tmax_gt90.missing_value = -32768\n",
    "        \n",
    "        temp_tmax_comfort = ncfile.createVariable('tmax_comfort',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmax_comfort.units = 'days' # number of days\n",
    "        temp_tmax_comfort.standard_name = 'estimated number of \"comfortable\" days where maximum daily temperature is '+ \\\n",
    "                                          'greater than 70 F and less than 80 F across 20-year window'\n",
    "        temp_tmax_comfort.missing_value = -32768\n",
    "        \n",
    "        temp_tmin_gt70 = ncfile.createVariable('tmin_gt70',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmin_gt70.units = 'days' # number of days\n",
    "        temp_tmin_gt70.standard_name = 'estimated number of days where minimum daily temperature equals or exceeds ' + \\\n",
    "                                       '70 F across 20-year window'\n",
    "        temp_tmin_gt70.missing_value = -32768\n",
    "        \n",
    "        temp_tmin_lt32 = ncfile.createVariable('tmin_lt32',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmin_lt32.units = 'days' # number of days\n",
    "        temp_tmin_lt32.standard_name = 'estimated number of days where minimum daily temperature equals or is less than '+ \\\n",
    "                                       '32 F across 20-year window'\n",
    "        temp_tmin_lt32.missing_value = -32768\n",
    "        \n",
    "        temp_tmin_lt0 = ncfile.createVariable('tmin_lt0',np.float32,('time','lat','lon')) # note: unlimited dimension is leftmost\n",
    "        temp_tmin_lt0.units = 'days' # number of days\n",
    "        temp_tmin_lt0.standard_name = 'estimated number of days where minimum daily temperature equals or is less than' + \\\n",
    "                                       '0 F across 20-year window'\n",
    "        temp_tmin_lt0.missing_value = -32768\n",
    "        \n",
    "\n",
    "        # Write latitudes, longitudes\n",
    "        # Note: the \":\" is necessary in these \"write\" statements\n",
    "        n_lat = n.variables['lat']\n",
    "        n_lon = n.variables['lon']\n",
    "        lat[:] = n_lat[:]\n",
    "        lon[:] = n_lon[:]\n",
    "\n",
    "        # write temp variables\n",
    "        temp_tmax_mean[:,:,:] = tmax_months_mean\n",
    "        temp_tmax_std[:,:,:] = tmax_months_std\n",
    "        temp_tmin_mean[:,:,:] = tmin_months_mean\n",
    "        temp_tmin_std[:,:,:] = tmin_months_std\n",
    "        temp_tmax_gt100[:,:,:] = tmax_months_gt100\n",
    "        temp_tmax_gt95[:,:,:] = tmax_months_gt95\n",
    "        temp_tmax_gt90[:,:,:] = tmax_months_gt90\n",
    "        temp_tmax_comfort[:,:,:] = tmax_months_comfort\n",
    "        temp_tmin_gt70[:,:,:] = tmin_months_gt70\n",
    "        temp_tmin_lt32[:,:,:] = tmin_months_lt32\n",
    "        temp_tmin_lt0[:,:,:] = tmin_months_lt0\n",
    "\n",
    "        # write time\n",
    "        yystart = y\n",
    "        nyears = 1\n",
    "        ntime = 12\n",
    "        print(\"Writing time. Year = \"+str(y))\n",
    "        datesout = [datetime.datetime(y,mm,15,0) for yy in range(yystart,yystart+nyears) for mm in range(1,13)]\n",
    "        time[:] = date2num(datesout,timeunits)\n",
    "\n",
    "        # close file\n",
    "        ncfile.close()\n",
    "\n",
    "        # report progress\n",
    "        print(\"Wrote: \"+filename)\n",
    "        \n",
    "        \n",
    "        # clear some memory, invoke Python garbage collector\n",
    "        del tmax_all_array\n",
    "        del tmax_means_array\n",
    "        del tmin_all_array\n",
    "        del tmin_means_array\n",
    "        del tmax_months_mean\n",
    "        del tmin_months_mean\n",
    "        del tmax_months_std\n",
    "        del tmin_months_std\n",
    "        del tmax_months_gt100\n",
    "        del tmax_months_gt95\n",
    "        del tmax_months_gt90\n",
    "        del tmax_months_comfort\n",
    "        del tmin_months_gt70\n",
    "        del tmin_months_lt32\n",
    "        del tmin_months_lt0\n",
    "        del tmax_mean\n",
    "        del tmin_mean\n",
    "        gc.collect()\n",
    "        print(\"Cleared some memory up...restarting loop\")\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46a26fa7",
   "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}