Commit 1e7afc3f authored by CEVAER's avatar CEVAER
Browse files

Added status attributes and percents attributes. Generate netcdf file even when bad status

parent 72b00de2
......@@ -109,7 +109,8 @@ def compute_centers(sat_file, resolution, api_data=None):
sat_bb = shapely.wkt.loads(df.iloc[0]["bounding_box"])
land, land_non_prep = land_mask()
center_st = center_status(sat_bb, atcf_point, land)
proportion_valid_data_around_center(sat_bb, atcf_point, land_non_prep)
percent_outside, percent_inside_island, percent_non_usable = \
proportion_valid_data_around_center(sat_bb, atcf_point, land_non_prep)
atcf_rad34_seq = df.iloc[0]["atcf_rad34_seq"]
atcf_rad34_neq = df.iloc[0]["atcf_rad34_neq"]
......@@ -147,8 +148,8 @@ def compute_centers(sat_file, resolution, api_data=None):
status = track_eye_sar(lon, lat, lon_eye_fg, lat_eye_fg, windspd_hetero_msk, rad=max(2 * atcf_rmw, r34 / 2),
atcf_rmw=atcf_rmw)
if status is not True: # and status != 30:
return None, status
#if status is not True: # and status != 30:
# return None, status
## 5. Reapply heterogeneity mask outside the eye ##
......@@ -170,7 +171,9 @@ def compute_centers(sat_file, resolution, api_data=None):
windspd_eye_masked = windspd.copy()
return {"lon": lon, "lat": lat, "windspd_eye_masked": windspd_eye_masked, "cm_cyc": cm_cyc,
return {"center_status": center_st, "percent_outside": percent_outside,
"percent_inside_island": percent_inside_island, "percent_non_usable": percent_non_usable,
"lon": lon, "lat": lat, "windspd_eye_masked": windspd_eye_masked, "cm_cyc": cm_cyc,
"lons_eye_shape": lons_eye_shape, "lats_eye_shape": lats_eye_shape, "lon_eye": lon_eye,
"lat_eye": lat_eye, "lon_eye_old": lon_eye_old, "lat_eye_old": lat_eye_old,
"lon_eye_fg": lon_eye_fg, "lat_eye_fg": lat_eye_fg, "lon_eye_hwnd": lon_eye_hwnd,
......@@ -210,15 +213,43 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
if api_data is None:
api_data = api_req(sat_file)
eye_dict, status = compute_centers(sat_file, resolution, api_data)
sat_date = owi.getDateFromFname(os.path.basename(sat_file))
if track_gdf is None:
track_gdf = track_data(api_data.iloc[0]["sid"], sat_date)
if status is not True: # and status != 30:
logger.error(f"Failed to find correct tc center ! Status {status}")
save_netcdf(sat_file=sat_file,
path=netcdf_path,
center_status=eye_dict["center_status"],
percent_outside=eye_dict["percent_outside"],
percent_inside_island=eye_dict["percent_inside_island"],
percent_non_usable=eye_dict["percent_non_usable"],
lon_center=eye_dict["lon_eye"],
lat_center=eye_dict["lat_eye"],
lons_eye_shape=eye_dict["lons_eye_shape"],
lats_eye_shape=eye_dict["lats_eye_shape"],
lon_eye_lwind=eye_dict["lon_eye_old"],
lat_eye_lwind=eye_dict["lat_eye_old"],
lon_eye_fg=eye_dict["lon_eye_fg"],
lat_eye_fg=eye_dict["lat_eye_fg"],
lon_eye_hwind=eye_dict["lon_eye_hwnd"],
lat_eye_hwind=eye_dict["lat_eye_hwnd"],
i_hwind=eye_dict["i_hwind"],
j_hwind=eye_dict["j_hwind"],
msk_lwind=eye_dict["msk_lwind"],
msk_hwind=eye_dict["msk_hwind"],
status=status,
storm_name=api_data.iloc[0]["cyclone_name"],
atcf_vmax=api_data.iloc[0]["vmax (m/s)"],
sid=api_data.iloc[0]["sid"],
atcf_filename=os.path.basename(api_data.iloc[0]["track_file"]),
lon_algo=eye_dict["lon"],
lat_algo=eye_dict["lat"],
resolution=resolution)
return status
else:
logger.info("Saving plot...")
sat_date = owi.getDateFromFname(os.path.basename(sat_file))
if track_gdf is None:
track_gdf = track_data(api_data.iloc[0]["sid"], sat_date)
p = create_control_plot(sat_file, eye_dict["lon"], eye_dict["lat"],
eye_dict["lons_eye_shape"], eye_dict["lats_eye_shape"],
......@@ -244,16 +275,34 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
sat_mission=api_data.iloc[0]["mission"],
acq_time=pd.to_datetime(api_data.iloc[0]["acquisition_start_time"]))
save_netcdf(sat_file, netcdf_path, eye_dict["lon_eye"], eye_dict["lat_eye"],
eye_dict["lons_eye_shape"], eye_dict["lats_eye_shape"],
eye_dict["lon_eye_old"], eye_dict["lat_eye_old"],
eye_dict["lon_eye_fg"], eye_dict["lat_eye_fg"], eye_dict["lon_eye_hwnd"],
eye_dict["lat_eye_hwnd"], eye_dict["i_hwind"], eye_dict["j_hwind"], eye_dict["msk_lwind"],
eye_dict["msk_hwind"], api_data.iloc[0]["cyclone_name"],
api_data.iloc[0]["vmax (m/s)"], api_data.iloc[0]["sid"],
os.path.basename(api_data.iloc[0]["track_file"]),
eye_dict["lon"], eye_dict["lat"],
resolution)
save_netcdf(sat_file=sat_file,
path=netcdf_path,
center_status=eye_dict["center_status"],
percent_outside=eye_dict["percent_outside"],
percent_inside_island=eye_dict["percent_inside_island"],
percent_non_usable=eye_dict["percent_non_usable"],
lon_center=eye_dict["lon_eye"],
lat_center=eye_dict["lat_eye"],
lons_eye_shape=eye_dict["lons_eye_shape"],
lats_eye_shape=eye_dict["lats_eye_shape"],
lon_eye_lwind=eye_dict["lon_eye_old"],
lat_eye_lwind=eye_dict["lat_eye_old"],
lon_eye_fg=eye_dict["lon_eye_fg"],
lat_eye_fg=eye_dict["lat_eye_fg"],
lon_eye_hwind=eye_dict["lon_eye_hwnd"],
lat_eye_hwind=eye_dict["lat_eye_hwnd"],
i_hwind=eye_dict["i_hwind"],
j_hwind=eye_dict["j_hwind"],
msk_lwind=eye_dict["msk_lwind"],
msk_hwind=eye_dict["msk_hwind"],
status=status,
storm_name=api_data.iloc[0]["cyclone_name"],
atcf_vmax=api_data.iloc[0]["vmax (m/s)"],
sid=api_data.iloc[0]["sid"],
atcf_filename=os.path.basename(api_data.iloc[0]["track_file"]),
lon_algo=eye_dict["lon"],
lat_algo=eye_dict["lat"],
resolution=resolution)
#save_center_plot(path, sat_file, eye_dict["lon"], eye_dict["lat"], eye_dict["windspd_eye_masked"],
# eye_dict["cm_cyc"], eye_dict["lons_eye_shape"], eye_dict["lats_eye_shape"],
# eye_dict["lon_eye"], eye_dict["lat_eye"],
......
......@@ -16,10 +16,11 @@ def research_contour(mask, lon, lat):
return lon_lwind, lat_lwind
def save_netcdf(sat_file, path, lon_center, lat_center, lons_eye_shape, lats_eye_shape,
def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_island,
percent_non_usable, lon_center, lat_center, lons_eye_shape, lats_eye_shape,
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, storm_name,
atcf_vmax, sid, atcf_filename, lon, lat, resolution):
lat_eye_hwind, i_hwind, j_hwind, msk_lwind, msk_hwind, status, storm_name,
atcf_vmax, sid, atcf_filename, lon_algo, lat_algo, resolution):
sat_file_sw = owi.L2CtoLight(sat_file, resolution=resolution, mode="sw")
ds = xr.open_dataset(sat_file_sw)
......@@ -42,15 +43,24 @@ def save_netcdf(sat_file, path, lon_center, lat_center, lons_eye_shape, lats_eye
ds_polar_north.attrs["Storm ID"] = sid
ds_polar_north.attrs["Track source file"] = atcf_filename
ds_polar_north.attrs["Source satellite file"] = os.path.basename(sat_file_sw)
ds_polar_north.attrs["Track point status"] = ",".join([str(ct_st.value) for ct_st in center_status])
ds_polar_north.attrs["Percent outside"] = percent_outside
ds_polar_north.attrs["Percent inside is land"] = percent_inside_island
#ds_polar_north.attrs["Percent outside + inside is land"] = percent_inside_island
ds_polar_north.attrs["Status"] = status
#ds_polar_north.attrs["% of area of 100km around track point outside acquisition bounding box"] = percent_outside
#ds_polar_north.attrs["% of area of 100km around track point inside acquisition bounding and on land"] = percent_inside_island
#ds_polar_north.attrs["% of area of 100km around track point inside " \
# "acquisition bounding and on land or outside of bounding box"] = percent_inside_island
lon_lwind, lat_lwind = research_contour(msk_lwind, lon, lat)
lon_lwind, lat_lwind = research_contour(msk_lwind, lon_algo, lat_algo)
lwind_contour_coords = []
for c_i in range(0, len(lon_lwind)):
lwind_contour_coords.append((lon_lwind[c_i], lat_lwind[c_i]))
lwind_contour = LineString(lwind_contour_coords)
ds_polar_north.attrs["Low wind research area"] = lwind_contour.wkt
lon_hwind, lat_hwind = research_contour(msk_hwind, lon, lat)
lon_hwind, lat_hwind = research_contour(msk_hwind, lon_algo, lat_algo)
hwind_contour_coords = []
for c_i in range(0, len(lon_lwind)):
hwind_contour_coords.append((lon_hwind[c_i], lat_hwind[c_i]))
......
......@@ -1032,18 +1032,20 @@ def proportion_valid_data_around_center(sat_bb, atcf_pt, land):
circle_poly = transform(aeqd_to_wgs84, buffer)
bb_intersection = circle_poly.intersection(sat_bb)
percent_bb_intesection = 100 - (bb_intersection.area / circle_poly.area) * 100
logger.info(f"{percent_bb_intesection}% of area of 100km around atcf track is outside acquisition bounding box")
percent_bb_intersection = 100 - (bb_intersection.area / circle_poly.area) * 100
logger.info(f"{percent_bb_intersection}% of area of 100km around atcf track is outside acquisition bounding box")
intersect_land = bb_intersection.intersection(land)
percent_intersect_land = (intersect_land.area / circle_poly.area) * 100
logger.info(f"{percent_intersect_land}% of area of 100km around atcf track point that is inside acquisition "
f"bounding box is land.")
percent_sum_non_usable = percent_intersect_land + percent_bb_intesection
percent_sum_non_usable = percent_intersect_land + percent_bb_intersection
logger.info(f"{percent_sum_non_usable}% of area of 100km around track point is not usable (land or"
f" outside of acquisition).")
return percent_bb_intersection, percent_intersect_land, percent_sum_non_usable
......
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