Commit 351c9782 authored by CEVAER's avatar CEVAER
Browse files

Fixed wrong lon/lat grid. Use numpy meshgrid. 3km instead of 12km

parent 0bc30207
......@@ -9,7 +9,7 @@ mkdir -p "$reportDir"
mkdir -p "$reportDir/logs"
mkdir -p "$reportDir/status"
gen_plot_fix.py -p "$outDir" -f "$owiFile" -r 12 > "$reportDir/logs/$owibase.log" 2>&1
gen_plot_fix.py -p "$outDir" -f "$owiFile" -r 3 > "$reportDir/logs/$owibase.log" 2>&1
status=$?
echo $status > "$reportDir/status/$owibase.status"
......
......@@ -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
center_status, proportion_valid_data_around_center, land_mask, distance_to_coast
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
......@@ -45,7 +45,7 @@ def api_req(sat_file):
return df
def track_data(sid, around_date):
def track_data(sid, around_date, base_url="https://cyclobs.ifremer.fr/app/"):
"""
Query CyclObs track API to retrieve track data 3 days around a given datetime.
......@@ -61,7 +61,7 @@ def track_data(sid, around_date):
"""
min_date = (around_date - datetime.timedelta(days=3)).strftime("%Y-%m-%d")
max_date = (around_date + datetime.timedelta(days=3)).strftime("%Y-%m-%d")
api_req = f"https://cyclobs.ifremer.fr/app/api/track?include_cols=all&sid={sid}&freq=50&min_date={min_date}&max_date={max_date}"
api_req = f"{base_url}api/track?include_cols=all&sid={sid}&freq=50&min_date={min_date}&max_date={max_date}"
logger.info(f"Getting data from Cyclobs track API. Req : {str(api_req)}")
track = pd.read_csv(api_req)
......@@ -109,6 +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)
dist_to_coast = distance_to_coast(atcf_point, land_non_prep)
logger.info(f"Distance to coast of track point is {dist_to_coast} meters.")
percent_outside, percent_inside_island, percent_non_usable = \
proportion_valid_data_around_center(sat_bb, atcf_point, land_non_prep)
......@@ -173,6 +175,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,
"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,
......@@ -225,6 +228,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
percent_outside=eye_dict["percent_outside"],
percent_inside_island=eye_dict["percent_inside_island"],
percent_non_usable=eye_dict["percent_non_usable"],
distance_to_coast=eye_dict["dist_to_coast"],
lon_center=eye_dict["lon_eye"],
lat_center=eye_dict["lat_eye"],
lons_eye_shape=eye_dict["lons_eye_shape"],
......@@ -281,6 +285,7 @@ def find_center_plot_save(path, sat_file, resolution, path_fix, netcdf_path, api
percent_outside=eye_dict["percent_outside"],
percent_inside_island=eye_dict["percent_inside_island"],
percent_non_usable=eye_dict["percent_non_usable"],
distance_to_coast=eye_dict["dist_to_coast"],
lon_center=eye_dict["lon_eye"],
lat_center=eye_dict["lat_eye"],
lons_eye_shape=eye_dict["lons_eye_shape"],
......
......@@ -20,7 +20,7 @@ def research_contour(mask, lon, lat):
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,
percent_non_usable, distance_to_coast, 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):
......@@ -57,11 +57,13 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
ds_polar_north["track_point_flag"] = xr.DataArray(data=",".join([str(ct_st.value) for ct_st in center_status]),
attrs={"long_name": "Track point flags",
"flag_values": "inside, outside, land",
"flag_values": "inside, outside",
"flag_meanings": "inside: track point is inside "
"satellite acquisition. outside : track"
" point is outside acquisition. land: "
"track point is on land."})
" point is outside acquisition."})
ds_polar_north["track_point_flag"] = xr.DataArray(data=distance_to_coast,
attrs={"long_name": "Distance of track point to nearest coast",
"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"})
......
......@@ -13,6 +13,7 @@ import scipy.interpolate as interp
import shapely.geometry as shp
from polylabel import polylabel
import scipy.signal as signal
from pyproj import Geod
from scipy.spatial import ConvexHull
from scipy.stats import norm
import matplotlib as mpl
......@@ -25,7 +26,7 @@ from numba import jit
from shapely import wkt
from enum import Enum
import cartopy
from shapely.ops import unary_union, transform
from shapely.ops import unary_union, transform, nearest_points
from shapely.prepared import prep
import pyproj
......@@ -1001,6 +1002,17 @@ def land_mask():
return land, land_geom
def distance_to_coast(point, land):
pts = nearest_points(point, land)
geod = Geod(ellps='WGS84')
fas, baz, distance = geod.inv(point.x, point.y, pts[1].x, pts[1].y)
if land.contains(point):
distance = -distance
return distance
def center_status(sat_bb, atcf_pt, land):
status = []
if sat_bb.contains(atcf_pt):
......@@ -1008,8 +1020,8 @@ def center_status(sat_bb, atcf_pt, land):
else:
status.append(CenterStatus.OUTSIDE)
if land.contains(atcf_pt):
status.append(CenterStatus.LAND)
#if land.contains(atcf_pt):
# status.append(CenterStatus.LAND)
return status
......
......@@ -146,7 +146,8 @@ def degrees2meters(lon, lat):
def control_plot_from_shp(sat_file_ll, eye_shape, center_pt, fg_pt, hwind_pt, lwind_pt, hwind_research, lwind_research,
atcf_vmax, track_gdf, cyclone_name, backend):
atcf_vmax, track_gdf, cyclone_name, backend,
opts={"width": 2600, "height": 2200, "legend_position": 'right', "fontscale": 4}):
sat_date = owi.getDateFromFname(os.path.basename(sat_file_ll))
high_cmap = cm.get_cmap("high_wind_speed")
#sat_file = owi.L2CtoLight(sat_file, mode="ll_gd")
......@@ -211,11 +212,13 @@ def control_plot_from_shp(sat_file_ll, eye_shape, center_pt, fg_pt, hwind_pt, lw
title=f"{cyclone_name}, ATCF VMAX: {atcf_vmax} m/s. "
f"Acquisition at {sat_date.strftime('%Y-%m-%d %H:%M:%S')}",
show_grid=True, active_tools=['wheel_zoom'],
xlim=(min_lon_meter, max_lon_meter), ylim=(min_lat_meter, max_lat_meter))
xlim=(min_lon_meter, max_lon_meter),
ylim=(min_lat_meter, max_lat_meter),
legend_position="right")
#p = p.redim.range(Longitude=min_max_lon, Latitude=min_max_lat)
if not matplot:
p.opts(width=2600, height=2200, legend_position='right', fontscale=4)
p.opts(**opts)
else:
p.opts(aspect=1, fig_inches=10, fig_bounds=(0, 0, 1, 1), legend_position='right')
......
......@@ -39,6 +39,27 @@ def custom_meshgrid(theta, rad, data_shape):
return coord_theta, coord_rad
def lon_lat_orientation(polar_ds):
polar_ds.lon.values[:] = np.flip(polar_ds.lon.values, axis=0)
polar_ds.lat.values[:] = np.flip(polar_ds.lat.values, axis=0)
t_shp = polar_ds.theta.shape
# Get the index in theta array corresponding to the 90° angle. (as theta[idx_90] =~ 90)
idx_90 = int(t_shp[0] / 4)
# Getting azimuth to north for center pixel
azim_heading = polar_ds.heading_angle.sel(theta=0, rad=0)
# Getting the closest index to the azim_heading azimuth
idx_azim = np.searchsorted(polar_ds.theta.values, azim_heading, side="left")
# Roll the lon and lat arrays of 90° + azim_angle to match the correct orientation.
polar_ds.lon.values[:] = np.roll(polar_ds.lon.values, idx_90 + idx_azim, axis=0)
polar_ds.lat.values[:] = np.roll(polar_ds.lat.values, idx_90 + idx_azim, axis=0)
return polar_ds
def to_polar(ds, lon_center, lat_center):
data_shape = ds["wind_speed"].shape
max_radius = int(np.sqrt(data_shape[0] ** 2 + data_shape[1] ** 2) / 2)
......@@ -54,9 +75,7 @@ def to_polar(ds, lon_center, lat_center):
rads_meter = rads * pix_dist
#print(rads_meter.min(), rads_meter.max())
#coord_theta, coord_rad = np.meshgrid(theta, rads_meter)
coord_theta, coord_rad = custom_meshgrid(theta, rads_meter, data_shape)
coord_theta, coord_rad = np.meshgrid(theta, rads_meter, indexing="ij")
geod = Geod(ellps='WGS84')
lons_polar, lats_polar, baz = geod.fwd(ds.lon.values * 0. + lon_center, ds.lat.values * 0. + lat_center,
......@@ -81,6 +100,8 @@ def to_polar(ds, lon_center, lat_center):
"rad": (["rad"], rads_meter)
})
polar_ds = lon_lat_orientation(polar_ds)
return polar_ds
......
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