Commit 043c78d5 authored by CEVAER's avatar CEVAER
Browse files

Added vmax and rmax to fix

parent 16b563c1
......@@ -134,7 +134,7 @@ def compute_centers(sat_file, resolution, api_data=None):
# Computing extra info useful to estimate if it probable to find an eye center into that acquisition.
land, land_non_prep = land_mask()
center_st = center_status(sat_bb, atcf_point, land)
center_st = center_status(sat_bb, atcf_point)
# Get smallest distance between atcf point and nearest coast.
dist_to_coast = distance_to_polyg(atcf_point, land_non_prep)
# Get smallest distance between atcf point and acquisition bounding box.
......@@ -314,8 +314,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
resolution=resolution)
return status
else:
logger.info("Saving plot...")
logger.info("Saving...")
#p = create_control_plot(sat_file, eye_dict["lon"], eye_dict["lat"],
# eye_dict["lons_eye_shape"], eye_dict["lats_eye_shape"],
# eye_dict["lon_eye"], eye_dict["lat_eye"],
......@@ -327,7 +326,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
#
#save_plot(sat_file, path, p, prepend="control")
radii_quad, radii_nmi_quad = save_netcdf(sat_file=sat_file,
radii_quad, radii_nmi_quad, percent_valid_quad = save_netcdf(sat_file=sat_file,
path=netcdf_path,
center_status=eye_dict["center_status"],
percent_outside=eye_dict["percent_outside"],
......@@ -365,6 +364,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
lat_eye=eye_dict["lat_eye"],
radii_quad=radii_quad,
radii_nmi_quad=radii_nmi_quad,
percent_valid_quad=percent_valid_quad,
cyclone_name=api_data.iloc[0]["cyclone_name"],
sid=api_data.iloc[0]["sid"],
sub_basin=api_data.iloc[0]["sub_basin"],
......
......@@ -70,10 +70,10 @@ def radius_for_wind(polar_ds, wind_mean, wind_knots, idx_max):
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
# if percent_invalid > 75:
# return np.nan, np.nan
return int(wind_mean.isel(rad=idx).item()), int(wind_mean.isel(rad=idx)["rad"].item())
return int(wind_mean.isel(rad=idx).item()), int(wind_mean.isel(rad=idx)["rad"].item()), 100 - percent_invalid
def compute_radii(ds_polar_north):
......@@ -90,9 +90,11 @@ def compute_radii(ds_polar_north):
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)}
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 = {}
percent_valid_quad = {}
for quad, sl in quadrants.items():
# Select the data in the current quadrant.
polar_quad = ds_polar_north.sel(theta=sl)
......@@ -112,24 +114,33 @@ def compute_radii(ds_polar_north):
rmax = polar_quad.isel(rad=idx_max)["rad"].item()
rmax_nmi = rmax * 0.000539957
# Percent invalid for rmax
count_total = polar_quad.isel(rad=idx_max).wind_speed.values.shape[0]
wind_nan_count = np.sum(np.where(np.isnan(polar_quad.isel(rad=idx_max).wind_speed.values), 1, 0))
percent_valid = 100 - (wind_nan_count / count_total * 100)
# 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]
percent_valid_radii = [percent_valid, percent_valid]
for w in wind_vals:
if vmax_kts < w:
wind_r = -1
radius_w = -1
radius_nmi = -1
percent_valid = -1
else:
wind_r, radius_w = radius_for_wind(polar_quad, wind_mean, w, idx_max)
wind_r, radius_w, percent_valid = radius_for_wind(polar_quad, wind_mean, w, idx_max)
radius_nmi = radius_w * 0.000539957
radii.append(radius_w)
radii_nmi.append(radius_nmi)
percent_valid_radii.append(percent_valid)
radii_quad[quad] = radii
radii_nmi_quad[quad] = radii_nmi
percent_valid_quad[quad] = percent_valid_radii
return radii_quad, radii_nmi_quad
return radii_quad, radii_nmi_quad, percent_valid_quad
def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_island,
......@@ -162,7 +173,7 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
ds_polar_north["track_vmax"] = xr.DataArray(data=float(atcf_vmax),
attrs={"long_name": "Track VMAX", "units": "m/s"})
ds_polar_north["cyclone_category"] = xr.DataArray(data=spd2cat(atcf_vmax, unit="m/s"),
attrs={"long_name": "Cyclone category based on track VMAX"})
attrs={"long_name": "Cyclone category based on track VMAX"})
ds_polar_north["eye_shape"] = xr.DataArray(data=eye_shape.wkt, attrs={"long_name": "Eye shape (WKT)"})
ds_polar_north["eye_center"] = xr.DataArray(data=Point(lon_center, lat_center).wkt,
attrs={"long_name": "TC eye center"})
......@@ -181,10 +192,10 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
"satellite acquisition. outside : track"
" point is outside acquisition."})
ds_polar_north["distance_to_coast"] = xr.DataArray(data=distance_to_coast,
attrs={"long_name": "Distance of track point to nearest coast. "
"Negative value means track point is "
"inside land.",
"units": "meters"})
attrs={"long_name": "Distance of track point to nearest coast. "
"Negative value means track point is "
"inside land.",
"units": "meters"})
ds_polar_north["distance_to_acquisition"] = xr.DataArray(data=distance_to_acq,
attrs={"long_name": "Distance of track point to nearest"
" acquisition border.",
......@@ -201,10 +212,11 @@ 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)
radii_quad, radii_nmi_quad, percent_valid_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)}
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():
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.",
......@@ -220,16 +232,21 @@ 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=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",
"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=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",
"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"percent_valid_radii_{quad}"] = xr.DataArray(data=percent_valid_quad[quad], attrs={
"long_name": f"Percent of valid data used for calculation of "
f"Vmax, Rmax, R64, R50, R34 for {quad} quadrant(s)",
"units": "percent"
})
# ds_polar_north.attrs["TC center status"] = status
......@@ -262,4 +279,4 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
ds_polar_north.to_netcdf(path=outfile)
logger.info(f"Wrote netCDF product to {outfile}")
return radii_quad, radii_nmi_quad
return radii_quad, radii_nmi_quad, percent_valid_quad
......@@ -8,7 +8,33 @@ import numpy as np
logger = logging.getLogger(__name__)
def save_fix_analysis(path, sat_file, atcf_filename, lon_eye, lat_eye, radii_quad, radii_nmi_quad, cyclone_name,
def get_vmax(radii_nmi_quad, percent_valid_quad):
max_vmax = -1
max_vmax_quad = None
for quad, radii in radii_nmi_quad.items():
if radii[0] > max_vmax:
max_vmax = radii[0]
max_vmax_quad = quad
if max_vmax == -1 or percent_valid_quad[max_vmax_quad][0] < 75:
return radii_nmi_quad["all"][0], radii_nmi_quad["all"][1]
return radii_nmi_quad[max_vmax_quad][0], radii_nmi_quad[max_vmax_quad][1]
def remove_floating_part(df):
# This to remove the floating part (.0) of the value.
df["rad1"] = df["rad1"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
df["rad2"] = df["rad2"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
df["rad3"] = df["rad3"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
df["rad4"] = df["rad4"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
df["wind speed"] = df["wind speed"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
df["MRD"] = df["MRD"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
return df
def save_fix_analysis(path, sat_file, atcf_filename, lon_eye, lat_eye, radii_quad, radii_nmi_quad, percent_valid_quad,
cyclone_name,
sid, sat_mission,
sub_basin, acq_time):
atcf_outfile = os.path.join(path, atcf_filename)
......@@ -38,7 +64,7 @@ def save_fix_analysis(path, sat_file, atcf_filename, lon_eye, lat_eye, radii_qua
fix.loc[0, 'BASIN'] = basin.upper()
fix.loc[0, 'CY'] = str(stormid).zfill(2)
fix.loc[0, 'date'] = owi.getDateFromFname(os.path.basename(sat_file)).strftime("%Y%m%d%H%M")
fix.loc[0, "fix format"] = 70
fix.loc[0, "fix format"] = "70"
fix.loc[0, "fix type"] = mis_to_type[sat_mission]
fix.loc[0, "center/intensity"] = "C"
fix.loc[0, 'lat'] = lat_eye
......@@ -54,6 +80,10 @@ def save_fix_analysis(path, sat_file, atcf_filename, lon_eye, lat_eye, radii_qua
fix.loc[0, "observation sources"] = sat_mission
fix.loc[0, "comments"] = "Synthetic Aperture Radar 3KM Wind Speed Analysis"
vmax, rmax = get_vmax(radii_nmi_quad, percent_valid_quad)
fix.loc[0, "wind speed"] = str(int(round(vmax)))
fix.loc[0, "MRD"] = str(int(round(rmax)))
fix_radii = pd.DataFrame(columns=fix.columns)
# quadrants = ["all", "NE", "SE", "SW", "NW"]
wind_vals = [64, 50, 34]
......@@ -65,14 +95,14 @@ def save_fix_analysis(path, sat_file, atcf_filename, lon_eye, lat_eye, radii_qua
new_row["rad2"] = ""
new_row["rad3"] = ""
new_row["rad4"] = ""
if not np.isnan(radii_nmi_quad["NE"][i_wnd + 2]) and radii_nmi_quad["NE"][i_wnd + 2] != -1:
new_row["rad1"] = round(radii_nmi_quad["NE"][i_wnd + 2])
if not np.isnan(radii_nmi_quad["SE"][i_wnd + 2]) and radii_nmi_quad["SE"][i_wnd + 2] != -1:
new_row["rad2"] = round(radii_nmi_quad["SE"][i_wnd + 2])
if not np.isnan(radii_nmi_quad["SW"][i_wnd + 2]) and radii_nmi_quad["SW"][i_wnd + 2] != -1:
new_row["rad3"] = round(radii_nmi_quad["SW"][i_wnd + 2])
if not np.isnan(radii_nmi_quad["NW"][i_wnd + 2]) and radii_nmi_quad["NW"][i_wnd + 2] != -1:
new_row["rad4"] = round(radii_nmi_quad["NW"][i_wnd + 2])
if percent_valid_quad["NE"][i_wnd + 2] > 75 and radii_nmi_quad["NE"][i_wnd + 2] != -1:
new_row["rad1"] = str(int(round(radii_nmi_quad["NE"][i_wnd + 2])))
if percent_valid_quad["SE"][i_wnd + 2] > 75 and radii_nmi_quad["SE"][i_wnd + 2] != -1:
new_row["rad2"] = str(int(round(radii_nmi_quad["SE"][i_wnd + 2])))
if percent_valid_quad["SW"][i_wnd + 2] > 75 and radii_nmi_quad["SW"][i_wnd + 2] != -1:
new_row["rad3"] = str(int(round(radii_nmi_quad["SW"][i_wnd + 2])))
if percent_valid_quad["NW"][i_wnd + 2] > 75 and radii_nmi_quad["NW"][i_wnd + 2] != -1:
new_row["rad4"] = str(int(round(radii_nmi_quad["NW"][i_wnd + 2])))
fix_radii = fix_radii.append(new_row)
if os.path.isfile(atcf_outfile):
......@@ -82,10 +112,8 @@ def save_fix_analysis(path, sat_file, atcf_filename, lon_eye, lat_eye, radii_qua
existing_fix["date"] = existing_fix["date"].apply(lambda x: x.strftime("%Y%m%d%H%M"))
existing_fix["start time"] = existing_fix["start time"].apply(lambda x: x.strftime("%Y%m%d%H%M"))
existing_fix["end time"] = existing_fix["end time"].apply(lambda x: x.strftime("%Y%m%d%H%M"))
existing_fix["rad1"] = existing_fix["rad1"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
existing_fix["rad2"] = existing_fix["rad2"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
existing_fix["rad3"] = existing_fix["rad3"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
existing_fix["rad4"] = existing_fix["rad4"].apply(lambda x: str(int(x)) if not np.isnan(x) else "")
existing_fix = remove_floating_part(existing_fix)
to_save = existing_fix.append(fix_radii, ignore_index=True)
to_save["start time"] = to_save["start time"].astype("string")
......
......@@ -23,6 +23,9 @@ matplot = False
def contour_indices(arr, lon):
"""
Retrieve 2D indexes of first shape found in arr
"""
arr[arr.mask] = 0
shapes = next(features.shapes(arr))
......
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