Skip to content
Snippets Groups Projects
Commit 5e0e0deb authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Reorganize the script `plot_snapshot.py` so that we can plot several 

snapshots in the same figure: useful for visualizing overlaps. So move
lines plotting a snapshot inside a new function called
snapshot and add option `--dashed`.

Remove the option `-e` of `plot_snaphot.py` and add directory argument
instead. The idea is that we now consider the three shapefiles as a
whole, they should always be alone in a dedicated directory.

Move the three shapefiles of `Tests/Input/Region_4/2006_01_02` to a
dedicated directory.

In `plot_snapshot.py`, netCDF variables for wind are now expected to be
called ugos and vgos.

Add a test for overlap of two different snapshots.
parent 6784855a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 157 additions and 129 deletions
Analysis/plot_snapshot.py
+ 146
− 129
View file @ 5e0e0deb
......@@ -20,7 +20,6 @@ mn for a 180° by 180° window.
"""
import argparse
import shapefile
import numpy as np
import matplotlib.pyplot as plt
......@@ -28,141 +27,159 @@ import netCDF4
import math
import itertools
import cartopy.crs as ccrs
import sys
parser = argparse.ArgumentParser()
parser.add_argument("-v", "--velocity", help = "plot velocity field",
action = "store_true")
parser.add_argument("-s", "--scale", default = 20, type = float,
help = "scale of arrows for the velocity field")
parser.add_argument("-g", "--grid", help = "plot grid", action = "store_true")
parser.add_argument("-w", "--window", help = "choose a limited plot window",
action = "store_true")
parser.add_argument("-l", "--light", help = "lighter plot",
action = "store_true")
parser.add_argument("-e", "--end_filename", help = "add character string after "
"extremum, outermost contour and max_speed_contour",
default = "")
args = parser.parse_args()
reader_extr = shapefile.Reader("extremum" + args.end_filename)
reader_outer = shapefile.Reader("outermost_contour" + args.end_filename)
try:
reader_m_s = shapefile.Reader("max_speed_contour" + args.end_filename)
except shapefile.ShapefileException:
print("Shapefile max_speed_contour_1 not found. "
"Max-speed contours will not be plotted.")
m_s_iterShapes = itertools.repeat(None)
else:
m_s_iterShapes = reader_m_s.iterShapes()
with netCDF4.Dataset("h.nc") as f:
longitude = f.variables["longitude"][:]
latitude = f.variables["latitude"][:]
if args.window:
l = input("llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = ? ")\
.split(sep = ",")
llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = [float(v) for v in l]
if urcrnrlon - llcrnrlon > 360:
sys.exit("bad values of urcrnrlon and llcrnrlon")
longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
# (in [llcrnrlon, llcrnrlon + 2 pi[)
lon_mask = longitude <= urcrnrlon
lat_mask = np.logical_and(latitude >= llcrnrlat, latitude <= urcrnrlat)
else:
lon_mask = np.ones(len(longitude), dtype = bool)
lat_mask = np.ones(len(latitude), dtype = bool)
plt.figure()
src_crs = ccrs.PlateCarree()
projection = ccrs.PlateCarree()
##projection = ccrs.NorthPolarStereo()
ax = plt.axes(projection = projection)
if args.grid:
lon_2d, lat_2d = np.meshgrid(longitude[lon_mask], latitude[lat_mask])
ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
marker = "+", color = "gray", linestyle = "None")
for shape_rec_extr, shape_outer, shape_m_s \
in zip(reader_extr, reader_outer.iterShapes(), m_s_iterShapes):
points = shape_rec_extr.shape.points[0]
if args.window:
points[0] += math.ceil((llcrnrlon - points[0]) / 360) * 360
# (in [llcrnrlon, llcrnrlon + 2 pi[)
if not args.window or (points[0] <= urcrnrlon
and llcrnrlat <= points[1] <= urcrnrlat):
if shape_rec_extr.record.cyclone == 0:
# Anti-cyclone
color = "red"
else:
color = "blue"
from os import path
lines = ax.plot(points[0], points[1], markersize = 10, color = color,
fillstyle = "none", transform = src_crs)
def snaphot(directory, dashed = False):
"""Plots extrema, outermost contours and max-speed contours.
if shape_rec_extr.record.valid == 1:
if args.light:
lines[0].set_marker("+")
directory: containing the three shapefiles.
dashed: boolean."""
reader_extr = shapefile.Reader(path.join(directory, "extremum"))
reader_outer = shapefile.Reader(path.join(directory, "outermost_contour"))
try:
reader_m_s = shapefile.Reader(path.join(directory, "max_speed_contour"))
except shapefile.ShapefileException:
print("Shapefile max_speed_contour not found. "
"Max-speed contours will not be plotted.")
m_s_iterShapes = itertools.repeat(None)
else:
m_s_iterShapes = reader_m_s.iterShapes()
for shape_rec_extr, shape_outer, shape_m_s \
in zip(reader_extr, reader_outer.iterShapes(), m_s_iterShapes):
points = shape_rec_extr.shape.points[0]
if args.window:
points[0] += math.ceil((llcrnrlon - points[0]) / 360) * 360
# (in [llcrnrlon, llcrnrlon + 2 pi[)
if not args.window or (points[0] <= urcrnrlon
and llcrnrlat <= points[1] <= urcrnrlat):
if shape_rec_extr.record.cyclone == 0:
# Anti-cyclone
color = "red"
else:
lines[0].set_marker("o")
elif not args.light:
# Invalid outermost contour
lines[0].set_marker("s")
if not args.light:
ax.annotate(str(shape_rec_extr.record.eddy_index),
projection.transform_point(points[0], points[1],
src_crs),
xytext = (3, 3), textcoords = "offset points")
if shape_outer.shapeType != shapefile.NULL and not args.light \
or shape_rec_extr.record.valid:
points = np.array(shape_outer.points)
lines = ax.plot(points[:, 0], points[:, 1], color = color,
color = "blue"
lines = ax.plot(points[0], points[1], markersize = 10,
color = color, fillstyle = "none",
transform = src_crs)
if not args.light:
if shape_rec_extr.record.valid == 0:
# Invalid outermost contour
lines[0].set_marker("s")
lines[0].set_fillstyle("none")
elif shape_m_s == None or shape_m_s.shapeType == shapefile.NULL:
lines[0].set_marker("x")
if shape_rec_extr.record.valid == 1:
if args.light:
lines[0].set_marker("+")
else:
lines[0].set_marker("o")
elif not args.light:
# Invalid outermost contour
lines[0].set_marker("s")
if shape_m_s != None and shape_m_s.shapeType != shapefile.NULL:
points = np.array(shape_m_s.points)
if shape_rec_extr.record.cyclone == 0:
# Anti-cyclone
color = "magenta"
else:
color = "cyan"
if not args.light:
ax.annotate(str(shape_rec_extr.record.eddy_index),
projection.transform_point(points[0], points[1],
src_crs),
xytext = (3, 3), textcoords = "offset points")
if shape_outer.shapeType != shapefile.NULL and not args.light \
or shape_rec_extr.record.valid:
points = np.array(shape_outer.points)
lines = ax.plot(points[:, 0], points[:, 1], color = color,
transform = src_crs)
if not args.light: lines[0].set_marker("o")
if args.velocity:
with netCDF4.Dataset("uv.nc") as f:
quiver_return = ax.quiver(longitude[lon_mask], latitude[lat_mask],
f["u"][0, lat_mask][:, lon_mask],
f["v"][0, lat_mask][:, lon_mask],
scale = args.scale, scale_units = "width",
transform = src_crs)
plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
coordinates = "figure")
ax.gridlines(draw_labels = True)
ax.coastlines()
plt.show()
if not args.light:
if shape_rec_extr.record.valid == 0:
# Invalid outermost contour
lines[0].set_marker("s")
lines[0].set_fillstyle("none")
elif shape_m_s == None \
or shape_m_s.shapeType == shapefile.NULL:
lines[0].set_marker("x")
else:
lines[0].set_marker("o")
if dashed: lines[0].set_linestyle("dashed")
if shape_m_s != None and shape_m_s.shapeType != shapefile.NULL:
points = np.array(shape_m_s.points)
if shape_rec_extr.record.cyclone == 0:
# Anti-cyclone
color = "magenta"
else:
color = "cyan"
lines = ax.plot(points[:, 0], points[:, 1], color = color,
transform = src_crs)
if not args.light: lines[0].set_marker("o")
if dashed: lines[0].set_linestyle("dashed")
if __name__ == "__main__":
import argparse
import sys
parser = argparse.ArgumentParser()
parser.add_argument("-v", "--velocity", help = "plot velocity field",
action = "store_true")
parser.add_argument("-s", "--scale", default = 20, type = float,
help = "scale of arrows for the velocity field")
parser.add_argument("-g", "--grid", help = "plot grid",
action = "store_true")
parser.add_argument("-w", "--window", help = "choose a limited plot window",
action = "store_true")
parser.add_argument("-l", "--light", help = "lighter plot",
action = "store_true")
parser.add_argument("-d", "--dashed", action = "store_true",
help = "dashed linestyle, useful for a second snapshot")
parser.add_argument("directory", help = "containing the three shapefiles")
args = parser.parse_args()
with netCDF4.Dataset("h.nc") as f:
longitude = f.variables["longitude"][:]
latitude = f.variables["latitude"][:]
if args.window:
l = input("llcrnrlon, llcrnrlat, urcrnrlon, "
"urcrnrlat = ? ").split(sep = ",")
llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = [float(v) for v in l]
if urcrnrlon - llcrnrlon > 360:
sys.exit("bad values of urcrnrlon and llcrnrlon")
longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
# (in [llcrnrlon, llcrnrlon + 2 pi[)
lon_mask = longitude <= urcrnrlon
lat_mask = np.logical_and(latitude >= llcrnrlat, latitude <= urcrnrlat)
else:
lon_mask = np.ones(len(longitude), dtype = bool)
lat_mask = np.ones(len(latitude), dtype = bool)
plt.figure()
src_crs = ccrs.PlateCarree()
projection = ccrs.PlateCarree()
##projection = ccrs.NorthPolarStereo()
ax = plt.axes(projection = projection)
if args.grid:
lon_2d, lat_2d = np.meshgrid(longitude[lon_mask], latitude[lat_mask])
ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
marker = "+", color = "gray", linestyle = "None")
snaphot(args.directory, args.dashed)
if args.velocity:
with netCDF4.Dataset("uv.nc") as f:
quiver_return = ax.quiver(longitude[lon_mask], latitude[lat_mask],
f["ugos"][0, lat_mask][:, lon_mask],
f["vgos"][0, lat_mask][:, lon_mask],
scale = args.scale, scale_units = "width",
transform = src_crs)
plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
coordinates = "figure")
ax.gridlines(draw_labels = True)
ax.coastlines()
plt.show()
Tests/short_tests.json
+ 11
− 0
View file @ 5e0e0deb
......@@ -302,5 +302,16 @@
]
],
"stdin_filename": "$input_dir/successive_overlap_periodic_nml.txt"
},
{
"args" : "$compil_prod_dir/test_successive_overlap",
"title" : "Successive_overlap_different_snapshots",
"description": "Overlap of different snapshots.",
"required":
[
["$input_dir/Region_4/2006_01_01/Snapshot", "Snapshot_1"],
["$input_dir/Region_4/2006_01_02/Snapshot", "Snapshot_2"]
],
"stdin_filename": "$input_dir/successive_overlap_nml.txt"
}
]
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