Commit a90ea28b authored by CEVAER's avatar CEVAER
Browse files

manage cases where radii wind not reached. Fixed total point count

parent 3a03d676
......@@ -26,10 +26,14 @@ def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max):
wind_masked[:idx_max] = np.ma.masked
idx = np.nanargmin(np.abs(wind_masked - wind_ms))
count_total = polar_ds.isel(rad=idx).wind_speed.count(dim="theta")
logger.debug(f"Target wind m/s: {wind_ms}, found wind value : {wind_mean.isel(rad=idx).item()}")
count_total = polar_ds.isel(rad=idx).wind_speed.values.shape[0]
wind_nan_count = np.sum(np.where(np.isnan(polar_ds.isel(rad=idx).wind_speed.values), 1, 0))
percent_invalid = wind_nan_count / count_total * 100
logger.debug(f"Total count points at {polar_ds.isel(rad=idx)['rad'].item()}m : {count_total}. "
f"Number of invalid point at same distance : {wind_nan_count}")
if percent_invalid > 75:
return np.nan, np.nan
......@@ -102,8 +106,12 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
# Compute and add radii to xarray dataset
quadrants = {"all": slice(0, 360), "NE": slice(0, 90), "SE": slice(90, 180), "SW": slice(180, 270), "NW": slice(270, 360)}
for quad, sl in quadrants.items():
print(quad)
polar_quad = ds_polar_north.sel(theta=sl)
wind_mean = polar_quad.wind_speed.mean(dim="theta")
invalid_x, invalid_y = np.where(ds_polar_north["mask_flag"].values > 0)
wind_masked = polar_quad["wind_speed"]
wind_masked[invalid_x, invalid_y] = np.nan
wind_mean = wind_masked.mean(dim="theta")
vmax = wind_mean.max(dim="rad").item()
vmax_kts = vmax * 1.94384
idx_max = np.where(wind_mean == vmax)[0][0]
......@@ -114,9 +122,15 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
radii = []
radii_nmi = []
for w in wind_vals:
wind_r, radius_w = radius_for_wind(ds_polar_north, wind_mean, w, idx_max)
if vmax_kts < w:
wind_r = -1
radius_w = -1
radius_nmi = -1
else:
wind_r, radius_w = radius_for_wind(ds_polar_north, wind_mean, w, idx_max)
radius_nmi = radius_w * 0.000539957
radii.append(radius_w)
radii_nmi.append(radius_w * 0.000539957)
radii_nmi.append(radius_nmi)
ds_polar_north[f"radii_ms_{quad}"] = xr.DataArray(data=[vmax, rmax, radii[0], radii[1], radii[2]], attrs={
"long_name": f"Vmax, Rmax, R64, R50, R34 for {quad} quadrant(s) wind speed mean.",
......@@ -126,7 +140,10 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
f"Rmax : radius of maximum wind speed, "
f"R64 : radius where wind speed equals 34 knots, "
f"R50 : radius where wind speed equals 50 knots, "
f"R34 : radius where wind speed equals 64 knots."
f"R34 : radius where wind speed equals 64 knots. "
f"NaN means that >75% of the data is not available"
f"for the radius calculation. -1 means that a correct radius cannot be found for that "
f"wind speed."
})
ds_polar_north[f"radii_kts_{quad}"] = xr.DataArray(data=[vmax_kts, rmax_nmi, radii_nmi[0], radii_nmi[1], radii_nmi[2]],
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment