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.

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

Add mask to edge of domain in smoke preprocessing #556

Merged
merged 2 commits into from
Dec 2, 2024
Merged
Changes from all commits
Commits
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
34 changes: 30 additions & 4 deletions ush/interp_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ def date_range(current_day):
print(f'Searching for interpolated RAVE for {current_day}')

fcst_datetime = dt.datetime.strptime(current_day, "%Y%m%d%H")

start_datetime = fcst_datetime - dt.timedelta(days=1, hours=1)

fcst_dates = pd.date_range(start=start_datetime, periods=24, freq='H').strftime("%Y%m%d%H")
Expand Down Expand Up @@ -167,6 +168,30 @@ def generate_regrider(rave_avail_hours, srcfield, tgtfield, weightfile, inp_file

return(regridder, use_dummy_emiss)


#Masks the edges of the data array to remove artifacts due to interpolation method on the top and bottom boundaries.
def mask_edges(data, mask_width=1):
"""
param data: numpy array, the data to mask
param mask_width: int, the width of the mask at each edge
return: numpy array, the masked data
"""
original_shape = data.shape
if mask_width < 1:
return data # No masking if mask_width is less than 1

# Mask top and bottom rows
data[:mask_width, :] = np.nan
data[-mask_width:, :] = np.nan

# Mask left and right columns
data[:, :mask_width] = np.nan
data[:, -mask_width:] = np.nan
assert data.shape == original_shape, "Data shape altered during masking."

return(data)


#process RAVE available for interpolation
def interpolate_rave(RAVE, rave_avail, rave_avail_hours, use_dummy_emiss, vars_emis, regridder,
srcgrid, tgtgrid, rave_to_intp, intp_dir, src_latt, tgt_latt, tgt_lont, cols, rows):
Expand Down Expand Up @@ -203,16 +228,17 @@ def interpolate_rave(RAVE, rave_avail, rave_avail_hours, use_dummy_emiss, vars_e
src_QA = xr.where(ds_togrid['FRE'] > 1000, src_rate, 0.0)
srcfield.data[...] = src_QA[0, :, :]
tgtfield = regridder(srcfield, tgtfield)
masked_tgt_data = mask_edges(tgtfield.data, mask_width=1)

if svar == 'FRP_MEAN':
Store_by_Level(fout, 'frp_avg_hr', 'Mean Fire Radiative Power', 'MW', '3D', '0.f', '1.f')
tgt_rate = tgtfield.data
tgt_rate = masked_tgt_data
fout.variables['frp_avg_hr'][0, :, :] = tgt_rate
print('=============after regridding===========' + svar)
print(np.sum(tgt_rate))
print(np.nansum(tgt_rate))
elif svar == 'FRE':
Store_by_Level(fout, 'FRE', 'FRE', 'MJ', '3D', '0.f', '1.f')
tgt_rate = tgtfield.data
tgt_rate = masked_tgt_data
fout.variables['FRE'][0, :, :] = tgt_rate
except (ValueError, KeyError) as e:
print(f"Error processing variable {svar} in {rave_file_path}: {e}")
Expand All @@ -221,4 +247,4 @@ def interpolate_rave(RAVE, rave_avail, rave_avail_hours, use_dummy_emiss, vars_e
except (OSError, IOError, RuntimeError, FileNotFoundError, TypeError, IndexError, MemoryError) as e:
print(f"Error reading NetCDF file {rave_file_path}: {e}")
else:
print(f"File not found or dummy emissions required: {rave_file_path}")
print(f"File not found or dummy emissions required: {rave_file_path}")