Commit b67401d7 authored by CEVAER's avatar CEVAER
Browse files

Added comments

parent ba6a46fd
#!/bin/bash
# Generate the fix analysis .dat and .netcdf for all the SAR acquisitions.
# Usage #
# -p for the path under which to save plots (required)
# -r for the path under which to save reports (required)
......@@ -23,8 +25,10 @@ done
script_dir=$(dirname $(readlink -f "$0")) # ie $CONDA_PREFIX/bin
# Retrieve the current SAR path (the commit can change the path).
proc_dir=$(dirname $(curl -s http://www.ifremer.fr/cerweb/shoc/shoc_dailyupdate_names.html | grep 'http-equiv="refresh"' | sed -r 's/.*URL=(.*)".*/\1/' | xargs pathUrl.py --schema=dmz -p))
# Retrieve the SAR acquisition file list from the most recent commit.
file_list="$proc_dir/sarwingL2X/ConcatprocessOk.txt"
# If line_number (-n) is not specified, updating all data
......
......@@ -12,8 +12,9 @@ logger = logging.getLogger(__name__)
if __name__ == "__main__":
description = """
Starts eye center research from a SAR file (L2C format). Saves control plot under the given path (-p).
A fix containing the position found is saved under <path>/fix.
Starts eye center research from a SAR file (L2C format).
A .dat fix containing the position found is saved under <path>/fix.
A .nc fix containing the data reprojected on a polar grid is saved under <path>/product
"""
parser = argparse.ArgumentParser(description=description)
......@@ -35,10 +36,13 @@ if __name__ == "__main__":
logger.setLevel(logging.INFO)
logging.basicConfig(level=logging.INFO)
# Prepare the directories
fix_path = os.path.join(args.path, "fix")
os.makedirs(fix_path, exist_ok=True)
netcdf_path = os.path.join(args.path, "product")
os.makedirs(netcdf_path, exist_ok=True)
# Main process
status = find_center_plot_save(args.path, args.sat_file, args.resolution, fix_path, netcdf_path, file_logging=True,
debug=args.debug)
......
#!/bin/bash
# A wrapper to write stdout to log file and to write the exit status to file too.
outDir=$1
owiFile=$2
reportDir=$3
......
......@@ -73,20 +73,43 @@ def track_data(sid, around_date, base_url="https://cyclobs.ifremer.fr/app/"):
def compute_centers(sat_file, resolution, api_data=None):
"""
Main function that prepares data and searches the cyclone eye center.
This includes opening the data, removing outliers, applying an heterogeneity mask,
finding the a high-wind-speed area barycenter, a low-wind-speed area barycenter, then
based on that low wind speed area, finding the eye shape and finally the eye_center as the centroid of that
eye shape.
Parameters
----------
sat_file: str The SAR acquisition file path (L2C format) to process.
resolution: int The resolution to use to find the eye center. (can be 1, 3, 5, 10, 12, 25, 40, 50)
api_data: pandas.DataFrame Optional. A pandas dataframe containing cyclObs API results for the given sat_file.
If not given, the function will make the request itself.
Returns
-------
dict
A python dict containg various information about the found center and intermediate data used in the processing.
"""
cm_cyc = make_cmap(['k', 'darkslategrey', 'teal', 'darkturquoise',
'turquoise', 'paleturquoise', 'lightcyan', 'azure', 'ivory', 'gold', 'darkorange', 'orangered',
'firebrick', 'm', 'magenta'])
sat_file_ref = sat_file
# Getting the _ll_gd version of the acquisition.
if resolution == 1:
sat_file = owi.L2CtoLight(sat_file)
else:
sat_file = owi.L2CtoLight(sat_file, resolution=resolution)
sat_file_low = owi.L2CtoL2M(sat_file, resolution=25)
#sat_file_low = owi.L2CtoL2M(sat_file, resolution=25)
lon, lat, windspd, winddir, nrcs, nd, FB, footprint_wkt = open_SAR_image(sat_file)
lon_lr, lat_lr, windspd_lr, winddir_lr, nrcs_lr, nd_lr, FB_lr, footprint_wkt_lr = open_SAR_image(sat_file_low)
#lon_lr, lat_lr, windspd_lr, winddir_lr, nrcs_lr, nd_lr, FB_lr, footprint_wkt_lr = open_SAR_image(sat_file_low)
# API req for that sat_file if not given as input.
if api_data is None:
df = api_req(sat_file_ref)
else:
......@@ -96,8 +119,10 @@ def compute_centers(sat_file, resolution, api_data=None):
logger.error("Couldn't find corresponding data in cyclobs API")
return None, 100
# Create shapely Point using the cyclone track point colocated to that SAR acquisition.
atcf_point = shapely.wkt.loads(df.iloc[0]["track_point"])
footprint = shapely.wkt.loads(footprint_wkt)
# Shapely polygon of SAR acquisition bounding box
sat_bb = shapely.wkt.loads(footprint_wkt)
atcf_lon = atcf_point.x
atcf_lat = atcf_point.y
logger.info(f"Atcf lon/lat {atcf_lon}/{atcf_lat}")
......@@ -107,16 +132,19 @@ def compute_centers(sat_file, resolution, api_data=None):
logger.error("Atcf rmw is empty.")
return None, 110
sat_bb = shapely.wkt.loads(df.iloc[0]["bounding_box"])
# 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)
# Get smallest distance between atcf point and nearest coast.
dist_to_coast = distance_to_polyg(atcf_point, land_non_prep)
dist_to_acq = distance_to_polyg(atcf_point, footprint.boundary, neg_if_inside=False)
# Get smallest distance between atcf point and acquisition bounding box.
dist_to_acq = distance_to_polyg(atcf_point, sat_bb.boundary, neg_if_inside=False)
logger.info(f"Distance to coast of track point is {dist_to_coast} meters.")
logger.info(f"Distance between aquisition border and track point is {dist_to_acq} meters.")
percent_outside, percent_inside_island, percent_non_usable = \
proportion_valid_data_around_center(sat_bb, atcf_point, land_non_prep)
# Get atcf 34 radii and mean them
atcf_rad34_seq = df.iloc[0]["atcf_rad34_seq"]
atcf_rad34_neq = df.iloc[0]["atcf_rad34_neq"]
atcf_rad34_swq = df.iloc[0]["atcf_rad34_swq"]
......@@ -133,8 +161,9 @@ def compute_centers(sat_file, resolution, api_data=None):
## I guess this is supposed to exclude rain signatures from the process. To do that, a gradient and a sliding
# windows to compute diff with neighbours of each pixel is used.
[windspd, nd, nrcs], mask_flt = mask_outlier_points(lon, lat, windspd, [nd, nrcs])
windspd_lr.mask = mask_flt
#windspd_lr.mask = mask_flt
# Use atcf track point position as first guess
lon_eye_fg = atcf_lon
lat_eye_fg = atcf_lat
......@@ -153,9 +182,6 @@ 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
## 5. Reapply heterogeneity mask outside the eye ##
[windspd, nd, nrcs], mask_h = mask_hetero(lon, lat, windspd, [nd, nrcs], FB, lon_eye_old, lat_eye_old,
......@@ -186,12 +212,40 @@ def compute_centers(sat_file, resolution, api_data=None):
"lat_eye_hwnd": lat_eye_hwnd, "msk_hwind": msk_hwind, "msk_lwind": msk_lwind,
"i_hwind": i_hwind, "j_hwind": j_hwind, "i_lwind": i_lwind, "j_lwind": j_lwind}, status
# return lon, lat, windspd_eye_masked, cm_cyc, lons_eye_shape, lats_eye_shape, lon_eye, lat_eye, lon_eye_old,\
# lat_eye_old, lon_eye_fg, lat_eye_fg, is_ok
def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api_data=None, track_gdf=None, debug=False,
file_logging=True):
"""
Function that calls the eye center research procedure and that pass the result to save functions (to .dat and to .nc
files).
Parameters
----------
path: str Base path to which resulting data will be saved.
sat_file: str Path of satellite file to process.
resolution: int Resolution to use to search for the eye center (can be 1, 5, 10, 12, 25, 40, 50)
path_fix: str Path to which the fix .dat files will be saved
netcdf_path: str Path to which the .nc files will be saved.
api_data: pandas.DataFrame Optional. Cyclobs API data result for the given sat_file. If not given, the function will
itself call the API (so an access to https://cyclobs.ifremer.fr is required in that case).
track_gdf: geopandas.GeoDataFrame Optional. Cyclobs API data result for the track associated to the cyclone obseved
by sat_file acquisition. If not given, the function will itself call the API (so an access to
https://cyclobs.ifremer.fr is required in that case).
debug: bool True : activate DEBUG logging, else it uses INFO.
file_logging: bool If True, configures the logger to log data to file.
Returns
-------
int
Status code of the processing. Status are :
100 = Cyclobs API could not be reached.
10 = SAR acquisition contains too few valid pixels.
20 = Maximum number of iteration was reached when searching eye center.
30 = The found eye center is too far from the high wind speed area barycenter
1 = Unmanaged error
0 = OK
"""
# Remove previously defined loggers
for handler in logging.root.handlers[:]:
logging.root.removeHandler(handler)
......@@ -219,9 +273,12 @@ 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)
#sat_date = owi.getDateFromFname(os.path.basename(sat_file))
# Track_gdf was previously used when building plots after processing, but plot creation was moved outside so its
# no longer needed. Could be removed later on.
#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}")
......@@ -269,6 +326,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
# api_data.iloc[0]["vmax (m/s)"], track_gdf)
#
#save_plot(sat_file, path, p, prepend="control")
radii_quad, radii_nmi_quad = save_netcdf(sat_file=sat_file,
path=netcdf_path,
center_status=eye_dict["center_status"],
......
......@@ -139,7 +139,9 @@ def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
input can be scalar or arrays.
"""
# convert decimal degrees to radians
[lon1, lat1, lon2, lat2] = [np.radians(lon1), np.radians(lat1), np.radians(lon2), np.radians(lat2)]
......@@ -273,14 +275,6 @@ def proc_research_eye(lon, lat, var, lon_eye_ref, lat_eye_ref, rad_research):
# Get lon/lat of the barycenter of the low wind area.
lon_eye_lwnd = lon[i_eye_l, j_eye_l]
lat_eye_lwnd = lat[i_eye_l, j_eye_l]
#
# i_close,j_close = np.ma.where(haversine(lon,lat,lon_eye_lwnd,lat_eye_lwnd)<=rad_research)
# i_lwind_new,j_lwind_new = np.ma.where((var<np.ma.min(var[i_close,j_close])+0.05*(np.ma.max(var[i_close,j_close])-np.ma.min(var[i_close,j_close])))& \
# (haversine(lon,lat,lon_eye_lwnd,lat_eye_lwnd)<=rad_research))
# i_lwind_new,j_lwind_new = i_lwind_new[~np.isnan(tonan_outliers((i_lwind_new*j_lwind_new).astype(float),m=3))],\
# j_lwind_new[~np.isnan(tonan_outliers((i_lwind_new*j_lwind_new).astype(float),m=3))]
# i_eye_l = int(np.floor(np.ma.median(i_lwind_new)))
# j_eye_l = int(np.floor(np.ma.median(j_lwind_new)))
# Iteratively, compute low wind barycenters around the first one, and stock them in an array. #
# If a value is found 2 times in the array, it is kept as the final position #
......@@ -997,6 +991,15 @@ class CenterStatus(Enum):
def land_mask():
"""
Get polygon from cartopy containing worldwide land.
Returns
-------
tuple
A shapely prepared geometry of the polygone and the 'raw' polygon.
"""
land_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')
land_polygons = list(land_10m.geometries())
land_geom = unary_union(land_polygons)
......@@ -1008,6 +1011,21 @@ def land_mask():
def distance_to_polyg(point, polyg, neg_if_inside=True):
"""
Return smallest distance in meters between a shapely point and a polygon.
Parameters
----------
point shapely.geometry.Point
polyg shapely.geometry.Polygon
neg_if_inside boolean: If True, returns negative value if the point is inside the polygone. If false, the
distance will only be positive.
Returns
-------
float
Distance in meters
"""
# nearest points will return 0 if you search the nearest point between a point and a polygon, the point given being
# inside the polygon. To prevent that, you can convert the Polygon to a linestring using polygon.boundary.
pts = nearest_points(point, polyg)
......@@ -1021,20 +1039,51 @@ def distance_to_polyg(point, polyg, neg_if_inside=True):
return distance
def center_status(sat_bb, atcf_pt, land):
def center_status(sat_bb, atcf_pt):
"""
Check if atcf point is inside or outisde satellite bounding box.
Parameters
----------
polyg shapely.geometry.Polygon: Satellite bounding box
point shapely.geometry.Point: atcf point
Returns
-------
list
list containing CenterStatus flags
"""
status = []
if sat_bb.contains(atcf_pt):
status.append(CenterStatus.INSIDE)
else:
status.append(CenterStatus.OUTSIDE)
#if land.contains(atcf_pt):
# status.append(CenterStatus.LAND)
return status
def proportion_valid_data_around_center(sat_bb, atcf_pt, land):
"""
Compute 3 percentages that help to know if the eye center can be found.
% of area of 100km around atcf track is outside acquisition bounding box
% of area of 100km around atcf track point that is inside acquisition bounding box is land.
% of area of 100km around track point is not usable (land or outside of acquisition).
Parameters
----------
sat_bb: shapely.geometry.Polygon satellite bounding box
atcf_pt: shapely.geometry.Point atcf point
land: shapely.geometry.Polygon polygon of land (coastline)
Returns
-------
tuple
Contains 3 percentages (floats)
% of area of 100km around atcf track is outside acquisition bounding box
% of area of 100km around atcf track point that is inside acquisition bounding box is land.
% of area of 100km around track point is not usable (land or outside of acquisition).
"""
radius = 100000
local_azimuthal_projection = f"+proj=aeqd +R=6371000 +units=m +lat_0={atcf_pt.y} +lon_0={atcf_pt.x}"
wgs84_to_aeqd = partial(
......
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