Commit be3693c7 authored by CEVAER's avatar CEVAER
Browse files

Added radii and more variable attributes

parent 4de010ce
......@@ -19,6 +19,17 @@ def research_contour(mask, lon, lat):
return lon_lwind, lat_lwind
def radius_for_wind(wind_mean, wind_knots, idx_max):
wind_ms = wind_knots * 0.514444
wind_masked = np.ma.asarray(wind_mean.values)
wind_masked[:idx_max] = np.ma.masked
idx = np.nanargmin(np.abs(wind_masked - wind_ms))
print(wind_masked.shape, wind_masked[idx-10:idx+10])
return wind_mean.isel(rad=idx).item(), wind_mean.isel(rad=idx)["rad"].item()
def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_island,
percent_non_usable, distance_to_coast, distance_to_acq,
lon_center, lat_center, lons_eye_shape, lats_eye_shape,
......@@ -82,6 +93,48 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
"which is unusable for analysis, in percent."
})
# 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)
wind_mean = polar_quad.wind_speed.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:
wind_r, radius_w = radius_for_wind(wind_mean, w, idx_max)
radii.append(wind_r)
radii_nmi.append(wind_r * 0.000539957)
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.",
"units": "m/s meters meters meters meters",
"description": f"Computed vmax and radii for {quad} quadrant(s). 5 values are available : "
f"Vmax : maximum wind speed on the wind profile, "
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."
})
ds_polar_north[f"radii_kts_{quad}"] = xr.DataArray(data=[vmax_kts, rmax_nmi, radii_nmi[0], radii_nmi[1], radii_nmi[2]],
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",
"description": f"Computed vmax and radii for {quad} quadrant(s). 5 values are available : "
f"Vmax : maximum wind speed on the wind profile, "
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."
})
# ds_polar_north.attrs["TC center status"] = status
# ds_polar_north.attrs["Area within 100km around track point outside acquisition bounding box in percent"] = percent_outside
# ds_polar_north.attrs["Area of 100km around track point inside acquisition bounding and on land in percent"] = percent_inside_island
......
......@@ -119,18 +119,17 @@ def to_polar(ds, lon_center, lat_center):
# Keeping only valid indices
nan_ar[valid_indices] = var_polar[valid_indices]
nan_ar = var_orientation(nan_ar, idx_azim, theta.shape)
vars_polar[var_name] = (["theta", "rad"], nan_ar)
vars_polar[var_name] = (["theta", "rad"], nan_ar, ds[var_name].attrs)
vars_polar["lon"] = (["theta", "rad"], lons_polar)
vars_polar["lat"] = (["theta", "rad"], lats_polar)
vars_polar["lon"] = (["theta", "rad"], lons_polar, ds.lon.attrs)
vars_polar["lat"] = (["theta", "rad"], lats_polar, ds.lat.attrs)
polar_ds = xr.Dataset(data_vars=vars_polar, coords={
"theta": (["theta"], theta),
"rad": (["rad"], rads_meter)
})
#polar_ds = lon_lat_orientation(polar_ds)
return polar_ds
......
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