Commit f4483111 authored by CEVAER's avatar CEVAER
Browse files

Added comments

parent b67401d7
......@@ -12,6 +12,22 @@ logger = logging.getLogger(__name__)
def research_contour(mask, lon, lat):
Return lon and lat coordinates of the circular shape contour contained in mask. The three input arrays
must be of the same size.
mask: numpy.ndarray 2D array containing a shape
lon: numpy.ndarray 2D array containing lon coordinates
lat: numpy.ndarry 2D array containing lat coordinates
Tuple containing coordinates of shape contour. Longitude coordinates at first index and latitude coordinates
at second index.
idx0, idx1 = contour_indices(mask, lon)
lon_lwind = lon[idx0, idx1]
lat_lwind = lat[idx0, idx1]
......@@ -19,7 +35,25 @@ def research_contour(mask, lon, lat):
return lon_lwind, lat_lwind
def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max, da_valid, da_all):
def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max):
Finds inside a polar projected dataset the radius at which the mean of the wind speed equals wind_knots.
The index of the found radius will always be superior to idx_max.
If the number of valid points to compute the mean is below 75% then np.nan is returned instead of the radius
and its associated wind speed.
polar_ds: xarray.Dataset Dataset containing SAR data in polar projection.
wind_mean: xarray.DataArray wind_speed variable meaned along theta axis in m/s.
wind_knots: int or float of the wind for which to find the radius.
idx_max: int index above which the found radius index must be.
the wind_speed at the found radius and the corresponding radius.
wind_ms = wind_knots * 0.514444
wind_masked =
......@@ -43,25 +77,42 @@ def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max, da_valid, da_all):
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))
Finds radii for each quadrant (all, NE, SE, SW, NW) and for 64, 50 and 34 wind speed (in knots).
ds_polar_north: xarray.Dataset SAR dataset in polar projection.
tuple(dict(quadrant: list(rad64, rad50, rad34)), dict(quadrant: list(rad64, rad50, rad34)))
radii data in meters and nmi
# 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():
# Select the data in the current quadrant.
polar_quad = ds_polar_north.sel(theta=sl)
# Set data to nan where mask indicates other than valid data.
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
# Mean the wind_speed.
wind_mean = wind_masked.mean(dim="theta")
# Find the vmax and rmax
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
# Loops to find the radii and store them into a dict contaning lists.
wind_vals = [64, 50, 34]
radii = [vmax, rmax]
radii_nmi = [vmax_kts, rmax_nmi]
......@@ -71,7 +122,7 @@ def compute_radii(ds_polar_north):
radius_w = -1
radius_nmi = -1
wind_r, radius_w = radius_for_wind(polar_quad, wind_mean, w, idx_max, da_valid, da_all)
wind_r, radius_w = radius_for_wind(polar_quad, wind_mean, w, idx_max)
radius_nmi = radius_w * 0.000539957
......@@ -87,6 +138,10 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
lon_eye_lwind, lat_eye_lwind, lon_eye_fg, lat_eye_fg, lon_eye_hwind,
lat_eye_hwind, i_hwind, j_hwind, msk_lwind, msk_hwind, status, storm_name,
atcf_vmax, sid, atcf_filename, lon_algo, lat_algo, resolution):
Re-project the SAR data to polar, add variables containing various scalar data (processing intermadiate values,
radii, etc) and saves the xarray Dataset to .nc file.
sat_file_sw = owi.L2CtoLight(sat_file, resolution=resolution, mode="sw")
ds = xr.open_dataset(sat_file_sw)
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