From 3fa2ed7588dafd9f49a36181c86fe63edca4f248 Mon Sep 17 00:00:00 2001 From: Liam Wedell Date: Tue, 26 Nov 2024 17:20:39 -0600 Subject: [PATCH] Add mask to edge of domain in smoke preprocessing --- ush/interp_tools.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/ush/interp_tools.py b/ush/interp_tools.py index cc9394786..02c7e6e3d 100755 --- a/ush/interp_tools.py +++ b/ush/interp_tools.py @@ -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") @@ -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): @@ -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}") @@ -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}") \ No newline at end of file + print(f"File not found or dummy emissions required: {rave_file_path}")