Commit 7916e20c authored by CEVAER's avatar CEVAER
Browse files

Fixed issue when computing rads

parent 2429a72f
......@@ -19,7 +19,7 @@ def research_contour(mask, lon, lat):
return lon_lwind, lat_lwind
def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max):
def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max, da_valid, da_all):
wind_ms = wind_knots * 0.514444
wind_masked = np.ma.asarray(wind_mean.values)
......@@ -30,15 +30,55 @@ def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max):
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}. "
f"Percent invalid : {percent_invalid}")
if percent_invalid > 75:
return np.nan, np.nan
return wind_mean.isel(rad=idx).item(), wind_mean.isel(rad=idx)["rad"].item()
return int(wind_mean.isel(rad=idx).item()), int(wind_mean.isel(rad=idx)["rad"].item())
def compute_radii(ds_polar_north):
da_valid = (["theta", "rad"], np.where(np.isnan(ds_polar_north.wind_speed.values), 0, 1))
da_all = (["theta", "rad"], np.ones(ds_polar_north.wind_speed.values.shape))
# 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)}
radii_quad = {}
radii_nmi_quad = {}
for quad, sl in quadrants.items():
polar_quad = ds_polar_north.sel(theta=sl)
invalid_x, invalid_y = np.where(polar_quad["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]
rmax = polar_quad.isel(rad=idx_max)["rad"].item()
rmax_nmi = rmax * 0.000539957
wind_vals = [64, 50, 34]
radii = [vmax, rmax]
radii_nmi = [vmax_kts, rmax_nmi]
for w in wind_vals:
if vmax_kts < w:
wind_r = -1
radius_w = -1
radius_nmi = -1
else:
wind_r, radius_w = radius_for_wind(polar_quad, wind_mean, w, idx_max, da_valid, da_all)
radius_nmi = radius_w * 0.000539957
radii.append(radius_w)
radii_nmi.append(radius_nmi)
radii_quad[quad] = radii
radii_nmi_quad[quad] = radii_nmi
return radii_quad, radii_nmi_quad
def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_island,
......@@ -106,35 +146,12 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
"which is unusable for analysis, in percent."
})
radii_quad, radii_nmi_quad = compute_radii(ds_polar_north)
# 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():
polar_quad = ds_polar_north.sel(theta=sl)
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]
rmax = polar_quad.isel(rad=idx_max)["rad"].item()
rmax_nmi = rmax * 0.000539957
wind_vals = [64, 50, 34]
radii = []
radii_nmi = []
for w in wind_vals:
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_nmi)
ds_polar_north[f"radii_ms_{quad}"] = xr.DataArray(data=[vmax, rmax, radii[0], radii[1], radii[2]], attrs={
ds_polar_north[f"radii_ms_{quad}"] = xr.DataArray(data=radii_quad[quad], attrs={
"long_name": f"Vmax, Rmax, R64, R50, R34 for {quad} quadrant(s) wind speed mean.",
"units": "m/s meters meters meters meters",
"description": f"Computed vmax and radii for {quad} quadrant(s). 5 values are available : "
......@@ -148,7 +165,7 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
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]],
ds_polar_north[f"radii_kts_{quad}"] = xr.DataArray(data=radii_nmi_quad[quad],
attrs={
"long_name": f"Vmax, Rmax, R64, R50, R34 for {quad} quadrant(s) wind speed mean.",
"units": "knots n.mi n.mi n.mi n.mi",
......@@ -189,3 +206,5 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
logger.info(f"Writing netCDF product to {outfile}")
ds_polar_north.to_netcdf(path=outfile)
logger.info(f"Wrote netCDF product to {outfile}")
return radii_quad, radii_nmi_quad
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