From d6ce19d526f07b843cbe5c79cc4e292181e6679f Mon Sep 17 00:00:00 2001 From: akrherz Date: Fri, 5 Mar 2021 11:15:09 -0600 Subject: [PATCH] add dumper per requested format #80 --- scripts/gridorder2/dumper.py | 162 +++++++++++++++++++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 scripts/gridorder2/dumper.py diff --git a/scripts/gridorder2/dumper.py b/scripts/gridorder2/dumper.py new file mode 100644 index 00000000..814b9a7c --- /dev/null +++ b/scripts/gridorder2/dumper.py @@ -0,0 +1,162 @@ +"""Generate the requested output.""" + +from pyiem.dep import read_env +from pyiem.util import get_dbconn +from tqdm import tqdm +import pandas as pd +from pandas.io.sql import read_sql + +HUCS = ( + "102400090102 102300020307 102400020604 102801020801 071000040902 " + "070600040608 070802040604 070802090402" +).split() +XREF = { + 91: 1, + 92: 2, + 93: 3, + 94: 4, + 95: 5, + 96: 6, + 97: 7, + 98: 9, + 99: 8, + 100: 11, + 101: 10, +} + + +def get_flowpath_lengths(pgconn): + """Load some metadata.""" + return read_sql( + "SELECT scenario, huc_12, fpath, bulk_slope," + "ST_LENGTH(geom) as len from flowpaths " + "WHERE scenario >= 91 and scenario <= 101 and huc_12 in %s" + "ORDER by scenario ASC, huc_12 ASC", + pgconn, + params=(tuple(HUCS),), + index_col=None, + ) + + +def compute_yearly(df, bymax=False, individual=False): + """Compute yearly 2011-2020 totals.""" + envs = [] + for _, row in tqdm(df.iterrows(), total=len(df.index)): + envfn = ( + f"/i/{row['scenario']}/env/{row['huc_12'][:8]}/" + f"{row['huc_12'][8:]}/{row['huc_12']}_{row['fpath']}.env" + ) + if bymax: + envdf = ( + read_env(envfn) + .sort_values("precip", ascending=False) + .groupby("year") + .nth(0) + .reset_index() + ) + else: + envdf = read_env(envfn).groupby("year").sum().copy().reset_index() + envdf["scenario"] = row["scenario"] + envdf["huc_12"] = row["huc_12"] + envdf["flowpath"] = row["fpath"] + envdf["delivery"] = envdf["sed_del"] / row["len"] + envs.append(envdf) + # We now have dataframe with yearly flowpath totals + envdf = pd.concat(envs) + key = ["scenario", "huc_12", "year"] + if individual: + key = ["scenario", "huc_12", "flowpath", "year"] + if individual: + genvdf = envdf.set_index(key) + else: + genvdf = envdf.groupby(by=key).mean().copy() + if bymax: + # add back in a good date, maybe + genvdf["date"] = envdf.groupby(by=key).nth(0)["date"] + return genvdf + + +def main(): + """Do great things.""" + pgconn = get_dbconn("idep") + df = get_flowpath_lengths(pgconn) + # df221 = df[df['len'] > 22.1] + # avg_len_full = df.groupby(["scenario", "huc_12"]).mean() + # avg_len_221 = df221.groupby(["scenario", "huc_12"]).mean() + # flowpath_counts_full = df.groupby(["scenario", "huc_12"]).count() + # flowpath_counts_221 = df221.groupby(["scenario", "huc_12"]).count() + # For all hillslopes + result_full = compute_yearly(df, bymax=True, individual=True) + # Only those over > 22.1m + # result_221 = compute_yearly(df[df['len'] > 22.1], individual=True) + + # Result time + rows = [] + for huc_12 in tqdm(HUCS): + for scenario in XREF: + fdf = df[(df["huc_12"] == huc_12) & (df["scenario"] == scenario)] + for flowpath in fdf["fpath"].values: + for year in range(2011, 2021): + key = (scenario, huc_12, flowpath, year) + if key not in result_full.index: + print(key) + continue + rows.append( + { + "Flowpath ID": flowpath, + "HUC12": huc_12, + "Situation": XREF[scenario], + # "Year": year, + # "Date for daily max rainfall": result_full.at[ + # (scenario, huc_12, year), "date"], + "Date for daily max rainfall": result_full.at[ + key, "date" + ], + "Flowpath length(m)": float( + fdf[(fdf["fpath"] == flowpath)]["len"] + ), + "Daily Rainfall (mm)": result_full.at[ + key, "precip" + ], + "Daily Runoff (mm)": result_full.at[key, "runoff"], + "Daily Delivery (kg/ha)": result_full.at[ + key, "delivery" + ], + "Daily Detachment (kg/ha)": result_full.at[ + key, "av_det" + ], + "Slope Gradient(%)": float( + fdf[(fdf["fpath"] == flowpath)]["bulk_slope"] + ) + * 100.0, + # "Number of flowpath": flowpath_counts_full.at[ + # (scenario, huc_12), 'len'], + # "Avg Flowpath Length(m)": avg_len_full.at[ + # (scenario, huc_12), "len"], + # "Daily Rainfall (mm)": result_full.at[ + # (scenario, huc_12, year), "precip"], + # "Daily Runoff (mm)": result_full.at[ + # (scenario, huc_12, year), "runoff"], + # "Daily Delivery (kg/ha)": result_full.at[ + # (scenario, huc_12, year), "delivery"], + # "Daily Detachment (kg/ha)": result_full.at[ + # (scenario, huc_12, year), "av_det"], + # "Number of flowpath >22.1": flowpath_counts_221.at[ + # (scenario, huc_12), 'len'], + # "Avg Flowpath Length(m) >22.1": avg_len_221.at[ + # (scenario, huc_12), "len"], + # "Daily Runoff (mm) >22.1": result_221.at[ + # (scenario, huc_12, year), "runoff"], + # "Daily Delivery (kg/ha) >22.1": result_221.at[ + # (scenario, huc_12, year), "delivery"], + # "Daily Detachment (kg/ha) >22.1": result_221.at[ + # (scenario, huc_12, year), "av_det"], + } + ) + resdf = pd.DataFrame(rows) + with pd.ExcelWriter("daily.xlsx") as writer: + resdf.to_excel(writer, "Single Flowpath-yearly", index=False) + + +if __name__ == "__main__": + main()