Commit 2706a024 authored by CEVAER's avatar CEVAER
Browse files

Added function to plot from Netcdf shapes

parent fe425adb
......@@ -13,12 +13,13 @@ import geoviews.tile_sources as gts
from bokeh.resources import INLINE
import matplotlib.cm as cm
import colormap_ext
import math
matplot = False
if matplot:
hv.extension('matplotlib')
else:
hv.extension('bokeh')
#if matplot:
# hv.extension('matplotlib')
#else:
# hv.extension('bokeh')
def contour_indices(arr, lon):
......@@ -128,8 +129,97 @@ def create_control_plot(sat_file, lon, lat, lons_eye_shape, lats_eye_shape, lon_
return p
def save_plot(sat_file, path, plot, prepend=""):
def save_plot(sat_file, path, plot, prepend="", backend="bokeh"):
hv.extension(backend)
hv.save(plot, os.path.join(path, f'{prepend}_{os.path.basename(sat_file)}.html'), resources=INLINE)
if backend == "bokeh":
hv.save(plot, os.path.join(path, f'{prepend}_{os.path.basename(sat_file)}.png'), resources=INLINE, fmt="png")
elif backend == "matplotlib":
hv.save(plot, os.path.join(path, f'{prepend}_{os.path.basename(sat_file)}.png'), resources=INLINE, fmt="png")
def degrees2meters(lon, lat):
x = lon * 20037508.34 / 180
y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180)
y = y * 20037508.34 / 180
return x, y
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):
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")
ds = xr.open_dataset(sat_file_ll)
min_max_lon = (round(float(ds["lon"].min()), 3) - 3, round(float(ds["lon"].max()), 3) + 3)
min_max_lat = (round(float(ds["lat"].min()), 3) - 3, round(float(ds["lat"].max()), 3) + 3)
min_lon_meter, min_lat_meter = degrees2meters(min_max_lon[0], min_max_lat[0])
max_lon_meter, max_lat_meter = degrees2meters(min_max_lon[1], min_max_lat[1])
vdim = hv.Dimension("wind_speed", range=(0, 80))
gvds = gv.Dataset(ds["wind_speed"].squeeze(), vdims=[vdim])
img = gvds.to(gv.Image, ['lon', 'lat'])
if backend == "bokeh":
img.opts(cmap=high_cmap, colorbar=True, tools=["hover"])
elif backend == "matplotlib":
img.opts(cmap=high_cmap)
fg_p = gv.Points([fg_pt], label="First guess (from interpolated ATCF track)")\
.opts(size=10, color="green", line_color="black", line_width=2, tools=["hover"])
hwnd_p = gv.Points([hwind_pt], label="High wind speed barycenter")\
.opts(size=10, color="#EAE8FF", line_color="black", line_width=2, tools=["hover"])
found_p = gv.Points([lwind_pt], label="Low wind speed barycenter").opts(size=10, color="#724CF9",
line_color="black",
line_width=2,
tools=["hover"])
final_p = gv.Points([center_pt], label="Center of eye shape")\
.opts(size=10, color="black", marker="s", line_color="grey", line_width=2, tools=["hover"])
#hwnd_area = gv.Points(list(zip(lon[i_hwnd, j_hwnd], lat[i_hwnd, j_hwnd])), label="High wind speed area")\
# .opts(color="grey", alpha=0.2)
eye = gv.Path(eye_shape, label="Eye shape")\
.opts(line_width=4, line_color="black", show_legend=True)
high_contour = gv.Path(hwind_research, label="High wind speed research area")\
.opts(line_width=3, show_legend=True, line_color="grey")
low_contour = gv.Path(lwind_research, label="Low wind speed research area").opts(line_width=3,
show_legend=True,
line_color="#EAE8FF")
track_plot = gv.Points(track_gdf, vdims=["wind_speed (m/s)", "date"], label="Cyclone ATCF interpolated track").opts(
color="wind_speed (m/s)", cmap=high_cmap, size=10, tools=["hover"], clim=(0, 80), line_color="black",
line_width=1.5)
track_path = copy.copy(track_gdf)
# Converting geom Points to LineString
track_path["geometry"] = [
LineString([val, track_gdf.loc[idx + 1, "geometry"]]) if idx < len(track_gdf["geometry"]) - 1 else LineString([val, val])
for idx, val in track_gdf["geometry"].iteritems()]
# Removing last line because its geom is useless
track_path = track_path.iloc[0:-1]
track_path_plot = gv.Path(track_path, vdims=["wind_speed (m/s)", "date"],
label="Cyclone ATCF interpolated track").opts(
color="black", tools=["hover"], line_width=3, show_legend=True)
t_plot = hv.Overlay(track_path_plot * track_plot)
p = (gts.EsriImagery() * img * t_plot * high_contour * low_contour *
eye * fg_p * hwnd_p * found_p * final_p).opts(tools=["hover", "save"],
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))
#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=3)
else:
p.opts(aspect=1, fig_inches=10, fig_bounds=(0, 0, 1, 1), legend_position='right')
return p
def save_center_plot(path, sat_file, lon, lat, windspd_eye_masked, cm_cyc, lons_eye_shape, lats_eye_shape, lon_eye,
......
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