<jupyter_start><jupyter_text>This script compares age metrics among models that use uniform versus heterogeneous porosity.  The goal is to evaluate if and how complexity influences age metrics across scales.  The script reads in the model_grid.csv file for each model and uses that to tie each particle's initial location to select categories, such as the HUC, NLCD, coarse fraction, etc.  The model_grid.csv file was created via verion 2 of the general models / GRTD notebooks, tied to PFJ's repo called "GenMod" on the USGS GitLab site.  This notebook is used to generate figure 9 of the FWP comparisons paper.<jupyter_code>__author__ = '<NAME>'
%matplotlib notebook
import os, sys
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib
import matplotlib.cm as cm
from datetime import datetime
import gdal
from gdal import ogr, osr
import gen_mod_functions as gm
import flopy as fp
import pickle
from ipywidgets import interact, Dropdown, Text
from IPython.display import display
from sklearn import linear_model
from sklearn.metrics import mean_squared_error, r2_score
from scipy import stats
import csv

import fit_parametric_distributions

try:
    import rasterio
except:
    print('Install rasterio into your python environment. Or, skip plotting of geotiffs at the end of this notebook')

modifier = False
def ReturnEntry(sender):
    modifier.value = intext.value<jupyter_output><empty_output><jupyter_text>Specify user input, including list of models to compare and which axes to plot them, plus attributes to analyze.<jupyter_code>#simulate_list = ['FWP1L_zK', 'FWP5L_zK', 'FWP5L_hK']  # list of models (w/ matching directory name) to analyze
#compare_list = ('FWP5L_zK', 'FWP5L_hK')  # compare these models directly; both uni_n and het_n versions.  Use a tuple to preserve order
evaluate = ('FWP5L_hK', 'FWP5L_zK')  # which two models to compare with color flooding of model factors.  1st minus 2nd, so order matters!
abrevdict = {'FWP5L_zK':'5zK', 'FWP5L_hK':'5hK'}

# This variable facilitates scenario testing (stopping or not at weak sinks/sources; porosity values & configuration, etc.)
#scenario_name = None  # set to None if no additional text was added to the modpath file names.
#scenario_name = 'passthroughsnk'  # an optional text string added to MP file names. MUST match value used in NB 01aa!!!
#scenario_name = 'het_n'  
scenario_name = 'het_n_1Dquad' #  1Dquad is the interpolation scheme

# when comparing results of the scenario to the base case, which weighting option to use? 'flux' or 'volume'
base_case = 'flux'

# for plotting 1:1 graphs
#x1, y1 = 0, 1  # index of names in the yet-to-be-populated eptlist variable.  (0=base case; 1 = scenario)
surf_aq_lays = 3  # deepest layer of the surficial aquifer.
nrow, ncol = 930, 650  # easier to hardcode this than waste the time to read-in a model just to get DIS info.

# for labeling RTD plots
label_by_model_name = True  # if false, will include zone info in label.

category = 'HUC8' 
age_cutoff = 65
minptl = 500  # The minimum number of particles for EACH model within each HUC (eg: if FWP5L has 2000, 
               # but FWP1L has 800, none get plotted in the 1:1 plots.  Still included in the RTD plots.)
               # 500-1000 seems reasonable for 1:1 plots as it includes only HUCs with really refined RTDs;
               # however, 100 - 500 seems more reasonable if we want to visualize spatial patterns because it
               # allows more HUCs to be plotted.  100-500 is based on visually inspecting RTD curves.
                
# Columns in the model_grid.csv file to keep.  Purge all others.
mg_columns = ['node_num', 'row', 'col', 'HUC6', 'HUC8', 'HUC10', 'HUC12', 'ibound', 'gage_id', 'coarse_flag', 
              'qu_atlas', 'catchment', 'ssurgo', 'stream_order', 'surfmat']
                
purge_hucs = [40602, 4060200, 406020000, 40602000000]  # all hucs for Lake Michigan      

HUCshpdict = {'HUC8':'E:/HUCS/WBD_4n7/WBD_HUC8_UTMft_FWPdomain.shp', 
              'HUC10':'E:/HUCS/WBD_4n7/WBD_HUC10_UTMft_FWPdomain.shp', 
              'HUC12':'E:/HUCS/WBD_4n7/WBD_HUC12_UTMft_FWPdomain.shp'}

HUCtiffdict = {'HUC8':'E:/HUCS/WBD_4n7/HUC8_UTMft_FWP.tiff', 
              'HUC10':'E:/HUCS/WBD_4n7/HUC10_UTMft_FWP.tiff',
              'HUC12':'E:/HUCS/WBD_4n7/HUC12_UTMft_FWP.tiff'}

dischargerastdict = {'FWP1L_zK': '../TIFFs/FWP1L_zK_SFRgain.tif', 'FWP5L_zK':'../TIFFs/FWP5L_zK_SFRgain.tif', 'FWP5L_hK':
                     '../TIFFs/FWP5L_hK_SFRgain.tif'}
lossrastdict = {'FWP1L_zK': '../TIFFs/FWP1L_zK_SFRloss.tif', 'FWP5L_zK':'../TIFFs/FWP5L_zK_SFRloss.tif', 'FWP5L_hK':
                     '../TIFFs/FWP5L_hK_SFRloss.tif'}<jupyter_output><empty_output><jupyter_text>Prep the script for the models to be analyzed<jupyter_code>homes = ['../Models']
fig_dir = '../Figures'

if not os.path.exists(fig_dir):
    os.mkdir(fig_dir)

mfpth = '../executables/MODFLOW-NWT_1.0.9/bin/MODFLOW-NWT_64.exe'
mp_exe_name = '../executables/modpath.6_0/bin/mp6x64.exe' 

dir_list = []
modlist = []
i = 0
r = 0

path_dict = {}
dfdict = {}
totp = {}
for home in homes:
    if os.path.exists(home):
        for dirpath, dirnames, filenames in os.walk(home):
            for f in filenames:
                if os.path.splitext(f)[-1] == '.nam':
                    mod = os.path.splitext(f)[0]
                    i += 1
                    #if mod in compare_list:
                    if mod in evaluate:
                        modlist.append(mod)
                        dir_list.append(dirpath)
                        r += 1
                        path_dict[mod] = dirpath
                               
print('    {} models read'.format(i))
print('These {} models will be analyzed: {}'.format(r, modlist))
# generate list of nam files:
nam_list = []
for pth in dir_list:
    model = os.path.normpath(pth).split(os.sep)[2]
    nam_file = '{}.nam'.format(model)
    nam_list.append(nam_file)
# setup dictionaries of the MODFLOW units for proper labeling of figures.
lenunit = {0:'undefined units', 1:'feet', 2:'meters', 3:'centimeters'}
timeunit = {0:'undefined', 1:'second', 2:'minute', 3:'hour', 4:'day', 5:'year'}

# Create dictionary of multipliers for converting model time units to days
time_dict = dict()
time_dict[0] = 1.0 # undefined assumes days, so enter conversion to days
time_dict[1] = 24 * 60 * 60
time_dict[2] = 24 * 60
time_dict[3] = 24
time_dict[4] = 1.0
time_dict[5] = 1.0
# Read-in model info and check max/min nlay & create list of DIS objects.  
# Assumes all else is the same among models (hnoflow, hdry, etc)
print ('Reading model information')
nlay_min = 100
nlay_max = 0

dis_objs = []
for i, model in enumerate(nam_list):
    nam_file = model
    model_ws = dir_list[i]
    
    fpmg = fp.modflow.Modflow.load(nam_file, model_ws=model_ws, exe_name=mfpth, version='mfnwt', 
                                   load_only=['DIS', 'BAS6', 'UPW', 'OC'], check=False)

    dis = fpmg.get_package('DIS')
    dis_objs.append(dis)
    bas = fpmg.get_package('BAS6')
    upw = fpmg.get_package('UPW')
    oc = fpmg.get_package('OC')

    delr = dis.delr
    delc = dis.delc
    nlay = dis.nlay
    nrow = dis.nrow
    ncol = dis.ncol
    bot = dis.getbotm()
    top = dis.gettop()

    hnoflo = bas.hnoflo
    ibound = np.asarray(bas.ibound.get_value())
    hdry = upw.hdry
    
    if nlay > nlay_max:
        nlay_max = nlay
    if nlay < nlay_min:
        nlay_min = nlay
        
    print('  .. done reading model {}'.format(i+1))

print ('   ... all done') 

print('minimum layers in a model:  {}'.format(nlay_min))
print('maximum layers in a model:  {}'.format(nlay_max))<jupyter_output>Reading model information<jupyter_text>Process endpoint files and merge with model_grid.csv<jupyter_code>def purge(ep_data):
    pre_Quaternary = ep_data.loc[ep_data.rt>=2.6e6]
    pre_Cretaceous = ep_data.loc[ep_data.rt>=66e6]
    preCambrian = ep_data.loc[ep_data.rt>=541e6]
    pre_earth = ep_data.loc[ep_data.rt>=4.6e9]

    print('\nFor your information:')
    print('{} particles were simulated as being older than Earth!'.format(preCambrian.shape[0]))
    print('{} particles were simulated as being PreCambrian in age.'.format(preCambrian.shape[0]))
    print('{} particles were simulated as being Cretaceous in age or older.'.format(pre_Cretaceous.shape[0]))
    print('{} particles were simulated as being pre-Quaternary in age.'.format(pre_Quaternary.shape[0]))
    
    ep_data = ep_data.loc[ep_data.rt<4.6e9]
    print('Purged particles older than earth')
    return(ep_data)
# read-in the model_grid.csv file for each model.  Then create a dataframe from each csv & pickle file.

dfdict = {}
totp = {}

for i, model in enumerate(evaluate):
#for i, model in enumerate(compare_list):
    g = os.path.join(path_dict[model], 'model_grid.csv')
    z = os.path.join(path_dict[model], 'zone_df.csv')
    try:
        df = pd.read_csv(g)
        df.ibound.replace(0, np.nan, inplace=True)
        df = df[df.ibound.notnull()]
        df = df[mg_columns]  # keep just the desired fields
        # re-calculate 2D cell number b/c node_num is computed differently for model_grid.csv and what's in the ept file.
        df['cellnum2d'] = df.row * ncol + df.col
        
        zone_df = pd.read_csv(z, index_col=0)
        for group in zone_df:
            print('\nReading-in the EPT file for {} in {}'.format(group, model))
            if scenario_name == None:
                print("Oops, I didn't really plan for this.  Stopping")
                sys.exit()
            else:
                p1 = os.path.join(path_dict[model], '{}_{}_{}.mpend'.format(model, base_case, group))
                p2 = os.path.join(path_dict[model], '{}_{}_{}_{}.mpend'.format(model, base_case, group, scenario_name))
                mnz1 = '{}_{}'.format(model, group)
                mnz2 = '{}_{}_{}'.format(model, group, scenario_name)
                eptu1 = fit_parametric_distributions.read_endpoints(p1, dis, time_dict)
                eptu2 = fit_parametric_distributions.read_endpoints(p2, dis, time_dict)
                
            eptu1 = purge(eptu1)
            eptu2 = purge(eptu2)
            eptu1['cellnum2d'] = (eptu1['Initial Row']-1) * ncol + (eptu1['Initial Column'] -1)  # -1 to convert to 0-based
            eptu2['cellnum2d'] = (eptu2['Initial Row']-1) * ncol + (eptu2['Initial Column'] -1)  # -1 to convert to 0-based
            eptu1_mg = eptu1.join(df, on='cellnum2d', lsuffix='_ept', rsuffix='_mg')
            eptu2_mg = eptu2.join(df, on='cellnum2d', lsuffix='_ept', rsuffix='_mg')
            eptu1_mg = eptu1_mg[eptu1_mg['Initial Layer'] <= surf_aq_lays]  # ensure that we're only analyzing Glacial!
            eptu2_mg = eptu2_mg[eptu2_mg['Initial Layer'] <= surf_aq_lays]  # ensure that we're only analyzing Glacial!

            #dfdict[model] = eptu_mg
            dfdict[mnz1] = eptu1_mg
            dfdict[mnz2] = eptu2_mg
            totp[mnz1] = eptu1_mg.rt.count()
            totp[mnz2] = eptu2_mg.rt.count()
            
            #x1 = eptlist[x1]  # pull out the file name after supplying the index
            #y1 = eptlist[y1]
            
    except (AttributeError, ValueError, IOError, IndexError):
        print('ERROR. THIS CODE BLOCK DID NOT COMPLETE. TROUBLE-SHOOT AND TRY AGAIN')
        print('The error occured while working on this model: {}'.format(model))
        raise SystemExit()

eptlist = list(dfdict.keys())
print('....done')
print(eptlist)
print(modlist)  # can't control the order of this one
print(evaluate)<jupyter_output>['FWP5L_hK_zones[1]', 'FWP5L_hK_zones[1]_het_n_1Dquad', 'FWP5L_zK_zones[1]', 'FWP5L_zK_zones[1]_het_n_1Dquad']
['FWP5L_hK', 'FWP5L_zK']
('FWP5L_hK', 'FWP5L_zK')<jupyter_text>Figure 9 of FWP comparison paper.  Finalized on 3/20/2019<jupyter_code># 1:1 plot of Young Fraction for complex (X) versus simple (Y)  MF models for uniform and heterogeneous porosities for all HUC scales

# That is, we're testing to see if adding heterogeneous porosity increases or decreases the mis-match between simple and complexe MF models.
# If rho, std, & ME all improve or degrade, then that may indicate if age metrics from simple MF models can be improved by using texture to inform porosity.

x1, y1, x2, y2 = eptlist[0], eptlist[2], eptlist[1], eptlist[3]  # columns in eptlist. In this case (x=complexe MF; y=simple MF)
hucs = ['HUC8', 'HUC10', 'HUC12']
nhucs = len(hucs)
fig, [ax1, ax2, ax3] = plt.subplots(nrows=1, ncols=3, figsize=(11, 4))
axisdict = {0:ax1, 1:ax2, 2:ax3}

for k, each in enumerate(hucs):
    # make sure we get ALL unique HUC IDs
    u1 = dfdict[eptlist[0]][each].unique()
    u2 = dfdict[eptlist[2]][each].unique()
    u = np.append(u1, u2)
    uniques = np.unique(u)
    # remove any HUCs listed in purge_hucs
    for h in purge_hucs:
        ind = np.where(uniques==h)
        uniques = np.delete(uniques, ind)
    axis = axisdict[k]

    yfhucdict = {}
    yfvalues = []
    skiphuc = []
    for j, cat_val in enumerate(uniques):  # each HUC ID
        yfmoddict = {}
        for i, md in enumerate(eptlist):  # each of the 3 FWP models
            mn = eptlist[i]
            df = dfdict[md].loc[dfdict[md][each]==cat_val]
            if df.rt.count() >= minptl:
                youngdf = df.loc[df.rt < age_cutoff]
                yf = youngdf.rt.count() / df.rt.count()
                yfmoddict[mn] = yf
                yfvalues.append(yf)
            else:
                skiphuc.append(cat_val)
                break
        if cat_val not in skiphuc:
            yfhucdict[cat_val] = yfmoddict    

    ddd = pd.DataFrame(yfhucdict).T
    Prho_y1, Pp_y1 = stats.pearsonr(np.array(ddd[x1]), np.array(ddd[y1]))
    ME_y1 = np.mean(np.array(ddd[x1]) - np.array(ddd[y1]))
    absdiff_y1 = np.array(ddd[x1]) - np.array(ddd[y1])
    absdiff_y1 = np.where(absdiff_y1 <0, absdiff_y1 * -1, absdiff_y1)
    MAE_y1 = np.mean(absdiff_y1)
    RMSE_y1 = np.sqrt(np.mean((np.array(ddd[x1]) - np.array(ddd[y1]))**2))
    mYFzK1 = np.mean(np.array(ddd[y1]))
    mYFhK1 = np.mean(np.array(ddd[x1]))
    
    Prho_y2, Pp_y2 = stats.pearsonr(np.array(ddd[x2]), np.array(ddd[y2]))
    ME_y2 = np.mean(np.array(ddd[x2]) - np.array(ddd[y2]))
    absdiff_y2 = np.array(ddd[x2]) - np.array(ddd[y2])
    absdiff_y2 = np.where(absdiff_y2 <0, absdiff_y2 * -1, absdiff_y2)
    MAE_y2 = np.mean(absdiff_y2)
    RMSE_y2 = np.sqrt(np.mean((np.array(ddd[x2]) - np.array(ddd[y2]))**2))
    mYFzK2 = np.mean(np.array(ddd[y2]))
    mYFhK2 = np.mean(np.array(ddd[x2]))
    
    #plotting
    mini = min(min(ddd[x1]), min(ddd[x2]), min(ddd[y1]), min(ddd[y2]))
    maxi = max(max(ddd[x1]), max(ddd[x2]), max(ddd[y1]), max(ddd[y2]))
    yrange = maxi-mini
    if k == 0:
        ddd.plot(kind='scatter', x=x1, y=y1, marker='o', color='green', label= 'uniform porosity', ax=axis)
        ddd.plot(kind='scatter', x=x2, y=y2, marker='^', color='blue', label = 'hetero. porosity', ax=axis)
        axis.text(mini, mini+yrange, r'$r$ = {:3.2f}'.format(Prho_y1), fontsize=9, color='green')
        axis.text(mini, mini+0.95*yrange, r'$Mean YF 5lzK:$ = {:3.2f}'.format(mYFzK1), fontsize=9, color='green')
        axis.text(mini, mini+0.90*yrange, r'$Mean YF 5lhK:$ = {:3.2f}'.format(mYFhK1), fontsize=9, color='green')
        axis.text(mini, mini+0.85*yrange, r'$r$ = {:3.2f}'.format(Prho_y2), fontsize=9, color='blue')
        axis.text(mini, mini+0.80*yrange, r'$Mean YF 5lzK:$ = {:3.2f}'.format(mYFzK2), fontsize=9, color='blue')
        axis.text(mini, mini+0.75*yrange, r'$Mean YF 5lhK:$ = {:3.2f}'.format(mYFhK2), fontsize=9, color='blue')

        axis.legend(bbox_to_anchor=(0.5, 0.00), loc='lower center', bbox_transform=fig.transFigure, ncol=2)  # move the legend
        # below the plots.  https://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot
    else:
        ddd.plot(kind='scatter', x=x1, y=y1, marker='o', color='green', ax=axis)
        ddd.plot(kind='scatter', x=x2, y=y2, marker='^', color='blue', ax=axis)
        axis.text(mini, mini+yrange, r'$r$ = {:3.2f}'.format(Prho_y1), fontsize=9, color='green')
        axis.text(mini, mini+0.95*yrange, r'$Mean YF 5lzK:$ = {:3.2f}'.format(mYFzK1), fontsize=9, color='green')
        axis.text(mini, mini+0.90*yrange, r'$Mean YF 5lhK:$ = {:3.2f}'.format(mYFhK1), fontsize=9, color='green')
        axis.text(mini, mini+0.85*yrange, r'$r$ = {:3.2f}'.format(Prho_y2), fontsize=9, color='blue')
        axis.text(mini, mini+0.80*yrange, r'$Mean YF 5lzK:$ = {:3.2f}'.format(mYFzK2), fontsize=9, color='blue')
        axis.text(mini, mini+0.75*yrange, r'$Mean YF 5lhK:$ = {:3.2f}'.format(mYFhK2), fontsize=9, color='blue')
    
    axis.plot((mini, maxi), (mini, maxi), 'k--')
    axis.set_xlabel(''), axis.set_ylabel('')
    axis.set_title(each)
    
fig.text(0.5, 0.14, 'Fraction of young water (<65 yr) in the complex (5-layer-hK) models', ha='center')
fig.text(0.068, 0.57, 'Fraction of young water (<65 yr) \nin the simpler (5-layer-zK) models', va='center', ha='center', rotation='vertical')
plt.subplots_adjust(top = .9, bottom=0.25)  # make room for legend

dst = 'YFrac_121_hetn_hK2zK_3HUCS'
#dst = 'TEST.png'
dst_pth = os.path.join(fig_dir, dst)
plt.savefig(dst_pth)
#plt.close()
dst_pth<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13725263279712996, 'fraction_numerical': 0.03136211086679783, 'mean_word_length': 3.21639424249817, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 26, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '304268', 'license': 'CC0-1.0', 'n_tokens_mistral': 6902, 'n_tokens_neox': 6452, 'n_words': 1589, 'path': 'Notebooks/01cc_121comp_het_n_Fig9.ipynb', 'repo_name': 'pjuckem/GRTD', 'size': '17287'}