Commit 80dbad5a authored by CEVAER's avatar CEVAER
Browse files

Added distance_to_acq to netcdf

parent 80cc37a8
......@@ -3,7 +3,7 @@ import numpy as np
from shapely import wkt
from tc_center.functions import open_SAR_image, mask_outlier_points, mask_hetero, track_eye_sar, center_eye, make_cmap,\
center_status, proportion_valid_data_around_center, land_mask, distance_to_coast
center_status, proportion_valid_data_around_center, land_mask, distance_to_polyg
from tc_center.analysis_product import save_netcdf
from tc_center.plotting import create_control_plot, save_plot
from tc_center.fix import save_fix_analysis
......@@ -84,8 +84,8 @@ def compute_centers(sat_file, resolution, api_data=None):
sat_file_low = owi.L2CtoL2M(sat_file, resolution=25)
lon, lat, windspd, winddir, nrcs, nd, FB = open_SAR_image(sat_file)
lon_lr, lat_lr, windspd_lr, winddir_lr, nrcs_lr, nd_lr, FB_lr = open_SAR_image(sat_file_low)
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)
if api_data is None:
df = api_req(sat_file_ref)
......@@ -97,6 +97,7 @@ def compute_centers(sat_file, resolution, api_data=None):
return None, 100
atcf_point = shapely.wkt.loads(df.iloc[0]["track_point"])
footprint = 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}")
......@@ -109,8 +110,10 @@ 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)
dist_to_coast = distance_to_coast(atcf_point, land_non_prep)
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)
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)
......@@ -175,7 +178,7 @@ def compute_centers(sat_file, resolution, api_data=None):
return {"center_status": center_st, "percent_outside": percent_outside,
"percent_inside_island": percent_inside_island, "percent_non_usable": percent_non_usable,
"dist_to_coast": dist_to_coast,
"dist_to_coast": dist_to_coast, "dist_to_acq": dist_to_acq,
"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,
......@@ -229,6 +232,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
percent_inside_island=eye_dict["percent_inside_island"],
percent_non_usable=eye_dict["percent_non_usable"],
distance_to_coast=eye_dict["dist_to_coast"],
distance_to_acq=eye_dict["dist_to_acq"],
lon_center=eye_dict["lon_eye"],
lat_center=eye_dict["lat_eye"],
lons_eye_shape=eye_dict["lons_eye_shape"],
......@@ -286,6 +290,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
percent_inside_island=eye_dict["percent_inside_island"],
percent_non_usable=eye_dict["percent_non_usable"],
distance_to_coast=eye_dict["dist_to_coast"],
distance_to_acq=eye_dict["dist_to_acq"],
lon_center=eye_dict["lon_eye"],
lat_center=eye_dict["lat_eye"],
lons_eye_shape=eye_dict["lons_eye_shape"],
......
......@@ -20,7 +20,8 @@ def research_contour(mask, lon, lat):
def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_island,
percent_non_usable, distance_to_coast, lon_center, lat_center, lons_eye_shape, lats_eye_shape,
percent_non_usable, distance_to_coast, distance_to_acq,
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, status, storm_name,
atcf_vmax, sid, atcf_filename, lon_algo, lat_algo, resolution):
......@@ -61,8 +62,14 @@ 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",
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.",
"units": "meters"})
ds_polar_north["percent_outside"] = xr.DataArray(data=percent_outside,
attrs={"long_name": "Area within 100km around track point outside"
" acquisition bounding box, in percent"})
......
......@@ -101,6 +101,10 @@ def open_SAR_image(slice_SAR, extra_vars=[]):
windspd.mask[np.where(np.isnan(lon))] = True
lon[np.where(np.isnan(lon))] = np.nanmin(lon)
lat[np.where(np.isnan(lat))] = np.nanmin(lat)
# extract acquisition bounding box
footprint = getattr(nc, "footprint")
# extract additional variables if existing
vars_extra = []
if extra_vars != []:
......@@ -125,9 +129,9 @@ def open_SAR_image(slice_SAR, extra_vars=[]):
filt = FB
if vars_extra != []:
return lon, lat, windspd, winddir, nrcs, nd, filt, vars_extra
return lon, lat, windspd, winddir, nrcs, nd, filt, footprint, vars_extra
else:
return lon, lat, windspd, winddir, nrcs, nd, filt
return lon, lat, windspd, winddir, nrcs, nd, filt, footprint
def haversine(lon1, lat1, lon2, lat2):
......@@ -1002,13 +1006,18 @@ def land_mask():
return land, land_geom
def distance_to_coast(point, land):
pts = nearest_points(point, land)
def distance_to_polyg(point, polyg, neg_if_inside=True):
# 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)
geod = Geod(ellps='WGS84')
print()
print("dist", point.x, point.y, pts[1].x, pts[1].y, pts[0].x, pts[0].y)
fas, baz, distance = geod.inv(point.x, point.y, pts[1].x, pts[1].y)
if land.contains(point):
distance = -distance
if neg_if_inside:
if polyg.contains(point):
distance = -distance
return distance
......
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