Skip to content

Commit

Permalink
remove calls to NT spillover routines
Browse files Browse the repository at this point in the history
  • Loading branch information
Scott Stewart committed Nov 4, 2024
1 parent 3ddf88d commit 29bea40
Show file tree
Hide file tree
Showing 4 changed files with 3 additions and 184 deletions.
102 changes: 1 addition & 101 deletions seaice_ecdr/initial_daily_ecdr.py
Original file line number Diff line number Diff line change
Expand Up @@ -988,111 +988,11 @@ def compute_initial_daily_ecdr_dataset(
| ecdr_ide_ds["invalid_ice_mask"].data[0, :, :]
)

# If all possible ocean grid cells have NaN siconc,
# then this field has no valid data
# NOTE: invalid_ice_mask includes non-ocean grid cells

# TODO: Consider skipping to all-missing if all are missing
# is_potential_seaice = ~set_to_zero_sic
# is_all_missing = np.all(np.isnan(cdr_conc_raw[is_potential_seaice]))

cdr_conc = cdr_conc_raw.copy()
cdr_conc[set_to_zero_sic] = 0
cdr_conc_pre_spillover = cdr_conc.copy()

"""
# NOTE: This section was coded up so we could test whether or not we
# want to use it.
# NOTE: Setting cdr_conc_raw to max of NT/BT and then applying
# land spillover algorithm results in too much false siconc
# near land.
# So, we will calculate the BT spillover with NT2 and BT algs
# and use that as a mask for cdr_seaice_conc
# TODO: Some minimum thresholding already occurs for bt_conc
# for bt_conc (and nt_conc?). I don't think the
# land_spillover algorithm should have to rely on this.
# Apply NT2/BT to bt_conc
l90c = get_land90_conc_field(
hemisphere=hemisphere,
resolution=tb_data.resolution,
ancillary_source=ancillary_source,
)
adj123 = get_adj123_field(
hemisphere=hemisphere,
resolution=tb_data.resolution,
ancillary_source=ancillary_source,
)
# Apply the NT2 land spillover algorithm to the bt_conc field
bt_conc_filtered = bt_conc.copy()
bt_conc_filtered[set_to_zero_sic] = 0
bt_spillover_applied_nt2 = apply_nt2_land_spillover(
conc=bt_conc_filtered,
adj123=adj123.data,
l90c=l90c.data,
anchoring_siconc=0.0,
affect_dist3=False,
)
# Apply the BT land spillover algorithm to the bt_conc_NT2 field
bt_spillover_applied_nt2_bt = coastal_fix(
conc=bt_spillover_applied_nt2.copy(),
missing_flag_value=np.nan,
land_mask=non_ocean_mask.data,
minic=10,
fix_goddard_bt_error=True,
)
# Return NaN for missing and for land
bt_spillover_applied_nt2_bt[non_ocean_mask.data] = np.nan
nt_conc_spillovered = nt_conc.copy()
nt_conc_spillovered[adj123.data == 1] -= 60
nt_conc_spillovered[adj123.data == 2] -= 40
nt_conc_spillovered[adj123.data == 3] -= 20
nt_conc_spillovered[nt_conc_spillovered < 0] = 0
cdr_conc_raw_adj = calc_cdr_conc(
bt_conc=bt_spillover_applied_nt2_bt,
nt_conc=nt_conc_spillovered,
cdr_conc_min_fraction=0.1,
cdr_conc_max_fraction=1.0,
)
"""

if land_spillover_alg == "BT_NT" and n_valid_tb_grid_cells != 0:
# The BT_NT algorithm requires separate consideration of the
# boostrap and NASATeam fields
bt_asCDRv4_conc = compute_bt_asCDRv4_conc(
bt_conc,
ecdr_ide_ds["bt_weather_mask"].data[0, :, :],
ecdr_ide_ds["invalid_ice_mask"].data[0, :, :],
ecdr_ide_ds["non_ocean_mask"].data,
)
nt_asCDRv4_conc = compute_nt_asCDRv4_conc(
nt_conc,
weather_mask=ecdr_ide_ds["nt_weather_mask"].data[0, :, :],
invalid_ice_mask=ecdr_ide_ds["invalid_ice_mask"].data[0, :, :],
)
cdr_conc = land_spillover(
cdr_conc=cdr_conc,
hemisphere=hemisphere,
tb_data=tb_data,
algorithm=land_spillover_alg,
land_mask=non_ocean_mask.data,
ancillary_source=ancillary_source,
bt_conc=bt_asCDRv4_conc,
nt_conc=nt_asCDRv4_conc,
bt_wx=ecdr_ide_ds["bt_weather_mask"].data[0, :, :],
fix_goddard_bt_error=True,
)
is_ntwx_CDRv4 = (nt_conc > 0) & ecdr_ide_ds["nt_weather_mask"].data[0, :, :]
cdr_conc[is_ntwx_CDRv4] = 0

elif n_valid_tb_grid_cells != 0:
# Not doing BT_NT
if n_valid_tb_grid_cells != 0:
cdr_conc = land_spillover(
cdr_conc=cdr_conc,
hemisphere=hemisphere,
Expand Down
1 change: 0 additions & 1 deletion seaice_ecdr/intermediate_daily.py
Original file line number Diff line number Diff line change
Expand Up @@ -721,7 +721,6 @@ def create_standard_ecdr_for_dates(
@click.option(
"--land-spillover-alg",
required=True,
# type=click.Choice(("BT_NT", "NT2", "ILS")),
type=click.Choice(get_args(LAND_SPILL_ALGS)),
)
@click.option(
Expand Down
3 changes: 1 addition & 2 deletions seaice_ecdr/multiprocess_intermediate_daily.py
Original file line number Diff line number Diff line change
Expand Up @@ -178,9 +178,8 @@ def multiprocess_intermediate_daily(
@click.option(
"--land-spillover-alg",
required=False,
# type=click.Choice(["BT_NT", "NT2", "ILS"]),
type=click.Choice(get_args(LAND_SPILL_ALGS)),
default="BT_NT",
default="NT2_BT",
)
@click.option(
"--ancillary-source",
Expand Down
81 changes: 1 addition & 80 deletions seaice_ecdr/spillover.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
import numpy.typing as npt
from pm_icecon.bt.compute_bt_ic import coastal_fix
from pm_icecon.land_spillover import apply_nt2_land_spillover
from pm_icecon.nt.compute_nt_ic import apply_nt_spillover
from pm_tb_data._types import Hemisphere
from scipy.ndimage import binary_dilation, generate_binary_structure, shift

Expand All @@ -14,7 +13,6 @@
get_adj123_field,
get_land90_conc_field,
get_non_ocean_mask,
get_nt_minic,
get_nt_shoremap,
)
from seaice_ecdr.tb_data import (
Expand All @@ -23,7 +21,7 @@
from seaice_ecdr.util import get_ecdr_grid_shape

NT_MAPS_DIR = Path("/share/apps/G02202_V5/cdr_testdata/nt_datafiles/data36/maps")
LAND_SPILL_ALGS = Literal["BT_NT", "NT2", "ILS", "ILSb", "NT2_BT"]
LAND_SPILL_ALGS = Literal["NT2", "ILS", "ILSb", "NT2_BT"]


def convert_nonocean_to_shoremap(*, is_nonocean: npt.NDArray):
Expand Down Expand Up @@ -237,9 +235,6 @@ def land_spillover(
algorithm: LAND_SPILL_ALGS,
land_mask: npt.NDArray,
ancillary_source: ANCILLARY_SOURCES,
bt_conc=None, # only needed if the BT or NT spillover are used
nt_conc=None, # only needed if the BT or NT spillover are used
bt_wx=None, # only needed if the BT or NT spillover are used
fix_goddard_bt_error: bool = False, # By default, don't fix Goddard bug
) -> npt.NDArray:
"""Apply the land spillover technique to the CDR concentration field."""
Expand Down Expand Up @@ -306,80 +301,6 @@ def land_spillover(

spillover_applied = spillover_applied_nt2_bt

elif algorithm == "BT_NT":
# TODO: This algorithm is complicated by the NT algorithm
# and should not be supported. It remains here...
# because otherwise it is hard to understand how it
# was implemented.
non_ocean_mask = get_non_ocean_mask(
hemisphere=hemisphere,
resolution=tb_data.resolution,
ancillary_source=ancillary_source,
)

# Bootstrap alg
bt_conc[bt_wx] = 0
bt_conc[bt_conc == 110] = np.nan
spillover_applied_bt = coastal_fix(
conc=bt_conc,
missing_flag_value=np.nan,
land_mask=land_mask,
minic=10,
fix_goddard_bt_error=fix_goddard_bt_error,
)

# NT alg
# Apply the NT to the nt_conc field
# Only apply that to the cdr_conc field if nt_spilled > bt_conc
shoremap = get_nt_shoremap(
hemisphere=hemisphere,
resolution=tb_data.resolution,
ancillary_source=ancillary_source,
)

minic = get_nt_minic(
hemisphere=hemisphere,
resolution=tb_data.resolution,
ancillary_source=ancillary_source,
)
nt_conc_for_ntspillover = nt_conc.copy()
is_negative_ntconc = nt_conc < 0
nt_conc_for_ntspillover[is_negative_ntconc] = np.nan

spillover_applied_nt = apply_nt_spillover(
conc=nt_conc_for_ntspillover,
shoremap=shoremap,
minic=minic,
match_CDRv4_cdralgos=True,
)
nt_asCDRv4 = (10.0 * spillover_applied_nt).astype(np.int16)
nt_asCDRv4[is_negative_ntconc] = -10
nt_asCDRv4[land_mask] = -100

# This section is for debugging; outputs bt and bt after spillover
# bt_asCDRv4 = spillover_applied_bt.copy()
# bt_asCDRv4[np.isnan(bt_asCDRv4)] = 110
# bt_asCDRv4.tofile('bt_afterspill_v5.dat')
# nt_asCDRv4.tofile('nt_afterspill_v5.dat')

# Note: be careful of ndarray vs xarray here!
is_nt_spillover = (
(spillover_applied_nt != nt_conc) & (~non_ocean_mask.data) & (nt_conc > 0)
)
use_nt_spillover = (spillover_applied_nt > bt_conc) & (spillover_applied_bt > 0)

spillover_applied = spillover_applied_bt.copy()
spillover_applied[use_nt_spillover] = spillover_applied_nt[use_nt_spillover]

# Here, we remove ice that the NT algorithm removed -- with conc
# < 15% -- regardless of the BT conc value
is_ntspill_removed = (
is_nt_spillover
& (spillover_applied_nt >= 0)
& (spillover_applied_nt < 14.5)
)
spillover_applied[is_ntspill_removed.data] = 0

elif algorithm == "ILS":
# Prepare the ILS array using the adj123 field to init ils_arr
ils_arr = np.zeros(cdr_conc.shape, dtype=np.uint8)
Expand Down

0 comments on commit 29bea40

Please sign in to comment.