diff --git a/stglib/core/qaqc.py b/stglib/core/qaqc.py index 0b5521eb..06d9ef66 100644 --- a/stglib/core/qaqc.py +++ b/stglib/core/qaqc.py @@ -360,10 +360,26 @@ def trim_mask(ds, var): def trim_maxabs_diff(ds, var): if var + "_maxabs_diff" in ds.attrs: - val = ds.attrs[var + "_maxabs_diff"] - cond = np.abs(np.ediff1d(ds[var], to_begin=0)) > ds.attrs[var + "_maxabs_diff"] - affected = cond.sum() - ds[var][cond] = np.nan + if "sample" in ds[var].dims: + cond = ( + np.abs(ds[var].diff(dim="sample", label="upper")) + >= ds.attrs[var + "_maxabs_diff"] + ) + affected = cond.sum() + val = ds.attrs[var + "_maxabs_diff"] + + # pad bads mask for first element along dim so not set to all NaNs + padding = xr.zeros_like(ds[var].isel({"sample": 0})).astype(bool) + cond = xr.concat([padding, cond], dim="sample") + ds[var] = ds[var].where(~cond) + + else: + val = ds.attrs[var + "_maxabs_diff"] + cond = ( + np.abs(np.ediff1d(ds[var], to_begin=0)) > ds.attrs[var + "_maxabs_diff"] + ) + affected = cond.sum() + ds[var][cond] = np.nan notetxt = f"Values filled where data increases or decreases by more than {val} units in a single time step; {affected} values affected. " diff --git a/stglib/vec/cdf2nc.py b/stglib/vec/cdf2nc.py index 1db6c00a..37976f87 100644 --- a/stglib/vec/cdf2nc.py +++ b/stglib/vec/cdf2nc.py @@ -51,13 +51,18 @@ def cdf_to_nc(cdf_filename, atmpres=False): # Add EPIC and CMG attributes ds = aqdutils.ds_add_attrs(ds, inst_type="VEC") - ds = associate_z_coord(ds) + # ds = associate_z_coord(ds) + + ds = fill_snr(ds) + + ds = fill_cor(ds) for var in ds.data_vars: # need to do this or else a "coordinates" attribute with value of "burst" hangs around ds[var].encoding["coordinates"] = None ds = qaqc.trim_min(ds, var) ds = qaqc.trim_max(ds, var) + ds = qaqc.trim_maxabs_diff(ds, var) ds = qaqc.trim_min_diff(ds, var) ds = qaqc.trim_min_diff_pct(ds, var) ds = qaqc.trim_max_diff(ds, var) @@ -422,6 +427,58 @@ def combine_vars(ds): return ds +def fill_snr(ds): + """ + Fill velocity data with corresponding beam snr value threshold + """ + if "snr_threshold" in ds.attrs: + Ptxt = str(ds.attrs["snr_threshold"]) + + for v in [ + "u_1205", + "v_1206", + "w_1204", + ]: + if v in ds: + for bm in ds["beam"].values: + ds[v] = ds[v].where(ds.snr.sel(beam=bm) > ds.attrs["snr_threshold"]) + + histtext = "Filled velocity data using snr threshold of {} for corresponding beam(s).".format( + Ptxt + ) + ds = utils.insert_note(ds, v, histtext) + + ds = utils.insert_history(ds, histtext) + + return ds + + +def fill_cor(ds): + """ + Fill velocity data with corresponding beam correlation value threshold + """ + if "cor_threshold" in ds.attrs: + Ptxt = str(ds.attrs["cor_threshold"]) + + for v in [ + "u_1205", + "v_1206", + "w_1204", + ]: + if v in ds: + for bm in ds["beam"].values: + ds[v] = ds[v].where(ds.cor.sel(beam=bm) > ds.attrs["cor_threshold"]) + + histtext = "Filled velocity data using cor threshold of {} for corresponding beam(s).".format( + Ptxt + ) + ds = utils.insert_note(ds, v, histtext) + + ds = utils.insert_history(ds, histtext) + + return ds + + def time_encoding(ds): """ensure we don't set dtypes uint for CF compliance""" if "units" in ds["time"].encoding: @@ -445,9 +502,12 @@ def time_encoding(ds): def var_encoding(ds): - for v in ["burst", "orientation", "cor", "amp"]: - ds[v].encoding["dtype"] = "int32" - ds[v].encoding["_FillValue"] = -2147483648 + for var in ds.data_vars: + if ds[var].dtype == "float64": + ds[var].encoding["dtype"] = "float32" + for var in ["burst", "orientation", "amp"]: + ds[var].encoding["dtype"] = "int32" + ds[var].encoding["_FillValue"] = -2147483648 return ds