Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add figure preprocessing scripts to precip variability metrics #1069

Merged
merged 10 commits into from
Mar 22, 2024
Merged

Add figure preprocessing scripts to precip variability metrics #1069

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
b09b66a
add figure scripts
acordonez Mar 12, 2024
669f9af
run precommit
acordonez Mar 12, 2024
d6c816e
move files
acordonez Mar 13, 2024
035f654
move files
acordonez Mar 13, 2024
ec70943
update prdday_to func
acordonez Mar 19, 2024
49e1dda
Merge branch 'main' into 1064_ao_fig_scripts
acordonez Mar 22, 2024
2ea0664
fix index
acordonez Mar 22, 2024
4d1472e
fix indices
acordonez Mar 22, 2024
c9d991f
update readme
acordonez Mar 22, 2024
a61707a
add *
acordonez Mar 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
update prdday_to func
  • Loading branch information
@acordonez
acordonez committed Mar 19, 2024
commit ec709436b2c0a86e070666c34f9ec886a547b826
78 changes: 30 additions & 48 deletions pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area_freq_mean_regional.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,18 @@
import numpy as np


def prdday_to_frq3hridx(prdday, frequency):
frq3hr = 1.0 / (float(prdday) * 8.0)
idx = (np.abs(frequency - frq3hr)).argmin()
return int(idx)


def prdday_to_frq1didx(prdday, frequency):
frq24hr = 1.0 / (float(prdday))
idx = (np.abs(frequency - frq24hr)).argmin()
def prdday_to_frqidx(prdday, frequency, ntd):
"""
Find frequency index from input period
Input
- prdday: period (day)
- frequency: frequency
- ntd: number of time steps per day (daily: 1, 3-hourly: 8)
Output
- idx: frequency index
"""
frq = 1.0 / (float(prdday) * ntd)
idx = (np.abs(frequency - frq)).argmin()
return int(idx)


Expand All @@ -37,6 +40,7 @@ def prdday_to_frq1didx(prdday, frequency):
if hr == "day":
acordonez marked this conversation as resolved.
Show resolved Hide resolved
frqs_forced = ["semi-annual", "annual"]
frqs_unforced = ["synoptic", "sub-seasonal", "seasonal-annual", "interannual"]
ntd = 1
elif hr == "3hr":
frqs_forced = ["semi-diurnal", "diurnal", "semi-annual", "annual"]
frqs_unforced = [
Expand All @@ -46,6 +50,7 @@ def prdday_to_frq1didx(prdday, frequency):
"seasonal-annual",
"interannual",
]
ntd = 8

infile = open(fname, "rb")
psdm = pickle.load(infile)
Expand Down Expand Up @@ -74,59 +79,36 @@ def prdday_to_frq1didx(prdday, frequency):
for frq in frqs:
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@msahn In this loop, lines 74-130, I would like to double-check that all the time scales are being correctly sliced. If you could review this section I would greatly appreciate it.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For semi-annual and annual, I use np.amax instead of np.nanmean to select maximum power within their range as below. This is to handle CMIP models that have different periods for semi-annual and annual cycles because they use different calendars (e.g., 360-day, 365-day, and Gregorian).

pcmdi_metrics/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py

Lines 433 to 440 in 035f654

elif frq == "semi-annual": # 180day=<pr=<183day
idx2 = prdday_to_frqidx(180, frequency, ntd)
idx1 = prdday_to_frqidx(183, frequency, ntd)
amfm = np.amax(am[idx1 : idx2 + 1])
elif frq == "annual": # 360day=<pr=<366day
idx2 = prdday_to_frqidx(360, frequency, ntd)
idx1 = prdday_to_frqidx(366, frequency, ntd)
amfm = np.amax(am[idx1 : idx2 + 1])

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@msahn Thank you for pointing this out! I have another question about the sub-daily timescale indices. If I'm interpreting this correctly, are we getting the average of data at frequencies equal to or larger than the sub-daily frequency?

amfm = np.nanmean(am[idx1 + 1 :])

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the sub-daily timescale, we average frequencies larger than 1 day (pr<1day). The frequency of 1day is included in the synoptic timescale (1day=<pr<20day). This information is written as comments in the code.

pcmdi_metrics/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py

Line 441 in 035f654

elif frq == "sub-daily": # pr<1day

print(frq)
if frq == "semi-diurnal": # pr=0.5day
idx = prdday_to_frq1didx(0.5, frequency)
idx = prdday_to_frqidx(0.5, frequency, ntd)
amfm = am[idx]
elif frq == "diurnal": # pr=1day
idx = prdday_to_frq1didx(1, frequency)
idx = prdday_to_frqidx(1, frequency, ntd)
amfm = am[idx]
if frq == "semi-annual": # 180day=<pr=<183day
if hr == "day":
idx2 = prdday_to_frq1didx(180, frequency)
idx1 = prdday_to_frq1didx(183, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(180, frequency)
idx1 = prdday_to_frq3hridx(183, frequency)
amfm = np.nanmean(am[idx1 : idx2 + 1])
idx2 = prdday_to_frqidx(180, frequency, ntd)
idx1 = prdday_to_frqidx(183, frequency, ntd)
amfm = np.amax(am[idx1 : idx2 + 1])
elif frq == "annual": # 360day=<pr=<366day
if hr == "day":
idx2 = prdday_to_frq1didx(360, frequency)
idx1 = prdday_to_frq1didx(366, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(360, frequency)
idx1 = prdday_to_frq3hridx(366, frequency)
amfm = np.nanmean(am[idx1 : idx2 + 1])
idx2 = prdday_to_frqidx(360, frequency, ntd)
idx1 = prdday_to_frqidx(366, frequency, ntd)
amfm = np.amax(am[idx1 : idx2 + 1])
elif frq == "sub-daily": # pr<1day
idx1 = prdday_to_frq1didx(1, frequency)
idx1 = prdday_to_frqidx(1, frequency, ntd)
amfm = np.nanmean(am[idx1 + 1 :])
elif frq == "synoptic": # 1day=<pr<20day
if hr == "day":
idx2 = prdday_to_frq1didx(1, frequency)
idx1 = prdday_to_frq1didx(20, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(1, frequency)
idx1 = prdday_to_frq3hridx(20, frequency)
idx2 = prdday_to_frqidx(1, frequency, ntd)
idx1 = prdday_to_frqidx(20, frequency, ntd)
amfm = np.nanmean(am[idx1 + 1 : idx2 + 1])
elif frq == "sub-seasonal": # 20day=<pr<90day
if hr == "day":
idx2 = prdday_to_frq1didx(20, frequency)
idx1 = prdday_to_frq1didx(90, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(20, frequency)
idx1 = prdday_to_frq3hridx(90, frequency)
idx2 = prdday_to_frqidx(20, frequency, ntd)
idx1 = prdday_to_frqidx(90, frequency, ntd)
amfm = np.nanmean(am[idx1 + 1 : idx2 + 1])
elif frq == "seasonal-annual": # 90day=<pr<365day
if hr == "day":
idx2 = prdday_to_frq1didx(90, frequency)
idx1 = prdday_to_frq1didx(365, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(90, frequency)
idx1 = prdday_to_frq3hridx(365, frequency)
idx2 = prdday_to_frqidx(90, frequency, ntd)
idx1 = prdday_to_frqidx(365, frequency, ntd)
amfm = np.nanmean(am[idx1 + 1 : idx2 + 1])
elif frq == "interannual": # 365day=<pr
if hr == "day":
idx2 = prdday_to_frq1didx(365, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(365, frequency)
idx2 = prdday_to_frqidx(365, frequency, ntd)
amfm = np.nanmean(am[: idx2 + 1])

psdmfm[frc][mip][dat][var][dom][frq] = amfm.tolist()
Expand Down
78 changes: 30 additions & 48 deletions pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_freq_mean_regional.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,18 @@
import xcdat as xc


def prdday_to_frq3hridx(prdday, frequency):
frq3hr = 1.0 / (float(prdday) * 8.0)
idx = (np.abs(frequency - frq3hr)).argmin()
return int(idx)


def prdday_to_frq1didx(prdday, frequency):
frq24hr = 1.0 / (float(prdday))
idx = (np.abs(frequency - frq24hr)).argmin()
def prdday_to_frqidx(prdday, frequency, ntd):
"""
Find frequency index from input period
Input
- prdday: period (day)
- frequency: frequency
- ntd: number of time steps per day (daily: 1, 3-hourly: 8)
Output
- idx: frequency index
"""
frq = 1.0 / (float(prdday) * ntd)
idx = (np.abs(frequency - frq)).argmin()
return int(idx)


Expand All @@ -42,6 +45,7 @@ def prdday_to_frq1didx(prdday, frequency):
if hr == "day":
frqs_forced = ["semi-annual", "annual"]
frqs_unforced = ["synoptic", "sub-seasonal", "seasonal-annual", "interannual"]
ntd = 1
elif hr == "3hr":
frqs_forced = ["semi-diurnal", "diurnal", "semi-annual", "annual"]
frqs_unforced = [
Expand All @@ -51,6 +55,7 @@ def prdday_to_frq1didx(prdday, frequency):
"seasonal-annual",
"interannual",
]
ntd = 8

frcs = ["forced", "unforced"]
vars = ["power", "rednoise", "sig95"]
Expand Down Expand Up @@ -109,83 +114,60 @@ def prdday_to_frq1didx(prdday, frequency):
psfm[frc][mip][dat][var] = {}
for frq in frqs:
if frq == "semi-diurnal": # pr=0.5day
idx = prdday_to_frq3hridx(0.5, frequency)
idx = prdday_to_frqidx(0.5, frequency, ntd)
fm = power.isel({"frequency": idx}).data
elif frq == "diurnal": # pr=1day
idx = prdday_to_frq3hridx(1, frequency)
idx = prdday_to_frqidx(1, frequency, ntd)
fm = power.isel({"frequency": idx}).data
if frq == "semi-annual": # 180day=<pr=<183day
if hr == "day":
idx2 = prdday_to_frq1didx(180, frequency)
idx1 = prdday_to_frq1didx(183, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(180, frequency)
idx1 = prdday_to_frq3hridx(183, frequency)
idx2 = prdday_to_frqidx(180, frequency, ntd)
idx1 = prdday_to_frqidx(183, frequency, ntd)
fm = (
power.isel({"frequency": slice(idx1, idx2 + 1)})
.mean("frequency")
.max("frequency")
.data
)
elif frq == "annual": # 360day=<pr=<366day
if hr == "day":
idx2 = prdday_to_frq1didx(360, frequency)
idx1 = prdday_to_frq1didx(366, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(360, frequency)
idx1 = prdday_to_frq3hridx(366, frequency)
idx2 = prdday_to_frqidx(360, frequency, ntd)
idx1 = prdday_to_frqidx(366, frequency, ntd)
fm = (
power.isel({"frequency": slice(idx1, idx2 + 1)})
.mean("frequency")
.max("frequency")
.data
)
elif frq == "sub-daily": # pr<1day
idx1 = prdday_to_frq1didx(1, frequency)
idx1 = prdday_to_frqidx(1, frequency, ntd)
fm = (
power.isel({"frequency": slice(0, idx1 + 1)})
.mean("frequency")
.data
)
elif frq == "synoptic": # 1day=<pr<20day
if hr == "day":
idx2 = prdday_to_frq1didx(1, frequency)
idx1 = prdday_to_frq1didx(20, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(1, frequency)
idx1 = prdday_to_frq3hridx(20, frequency)
idx2 = prdday_to_frqidx(1, frequency, ntd)
idx1 = prdday_to_frqidx(20, frequency, ntd)
fm = (
power.isel({"frequency": slice(idx1 + 1, idx2 + 1)})
.mean("frequency")
.data
)
elif frq == "sub-seasonal": # 20day=<pr<90day
if hr == "day":
idx2 = prdday_to_frq1didx(20, frequency)
idx1 = prdday_to_frq1didx(90, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(20, frequency)
idx1 = prdday_to_frq3hridx(90, frequency)
idx2 = prdday_to_frqidx(20, frequency, ntd)
idx1 = prdday_to_frqidx(90, frequency, ntd)
fm = (
power.isel({"frequency": slice(idx1 + 1, idx2 + 1)})
.mean("frequency")
.data
)
elif frq == "seasonal-annual": # 90day=<pr<365day
if hr == "day":
idx2 = prdday_to_frq1didx(90, frequency)
idx1 = prdday_to_frq1didx(365, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq3hridx(90, frequency)
idx1 = prdday_to_frq3hridx(365, frequency)
idx2 = prdday_to_frqidx(90, frequency, ntd)
idx1 = prdday_to_frqidx(365, frequency, ntd)
fm = (
power.isel({"frequency": slice(idx1 + 1, idx2 + 1)})
.mean("frequency")
.data
)
elif frq == "interannual": # 365day=<pr
if hr == "day":
idx2 = prdday_to_frq1didx(365, frequency)
elif hr == "3hr":
idx2 = prdday_to_frq1didx(365, frequency)
idx2 = prdday_to_frqidx(365, frequency, ntd)
fm = (
power.isel({"frequency": slice(0, idx2 + 1)})
.mean("frequency")
Expand Down
Loading