From f88b44a4f239b1edcc650e3cb735ad5f2d5638a8 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Tue, 8 Jan 2019 16:45:43 +0100
Subject: [PATCH] Add default value of argument end_filename of function def
read in script distribution_function.py. Return dictionary instead of tuple:
easier to use.
In script plot_snapshot.py, move definition of readers to the
beginning so there is no risk of entering the window for nothing if we
forgot the end_filename option. Also, enter the four values defining
the window in a single input: easier to copy and paste.
---
Analysis/distribution_function.py | 43 ++++++++++++----------
Analysis/plot_snapshot.py | 45 ++++++++++++-----------
Documentation_texfol/documentation.tex | 51 ++++++++++++++++++++++++++
spherical_polygon_area.f | 6 +--
spherical_polyline_area.f | 8 ++--
5 files changed, 104 insertions(+), 49 deletions(-)
diff --git a/Analysis/distribution_function.py b/Analysis/distribution_function.py
index c1b7e5b3..68797f53 100755
--- a/Analysis/distribution_function.py
+++ b/Analysis/distribution_function.py
@@ -15,9 +15,10 @@ def plot_distr_funct(X, label = None, ax = None):
print("maximum value:", X[-1])
print("number of values:", nx, "\n")
-def read(end_filename):
- """Return speed, radius of outermost contour and radius of maximum speed
- contour from three dbf files, as Numpy arrays.
+def read(end_filename = ""):
+ """Read the three dbf files and return speed, radius and amplitude of
+ outermost contour, radius and amplitude of maximum speed contour,
+ as Numpy arrays.
Select valid speed values and valid outermost contour (with
sufficient area).
@@ -46,8 +47,9 @@ def read(end_filename):
rad_speed.append(rec_max.r_eq_area)
amp_speed.append(rec_extr.ssh - rec_max.ssh)
- return np.array(speed), np.array(rad_outer), np.array(rad_speed), \
- np.array(amp_outer), np.array(amp_speed)
+ return {"speed": np.array(speed), "rad_outer": np.array(rad_outer),
+ "rad_speed": np.array(rad_speed), "amp_outer": np.array(amp_outer),
+ "amp_speed": np.array(amp_speed)}
if __name__ == "__main__":
import argparse
@@ -58,62 +60,62 @@ if __name__ == "__main__":
"extremum, outermost_contour and max_speed_contour",
default = "")
args = parser.parse_args()
-
- speed, rad_outer, rad_speed, amp_outer, amp_speed = read(args.end_filename)
+ d = read(args.end_filename)
xlabel = "speed, in m s-1"
plt.figure()
- plt.hist(speed, bins = "auto")
+ plt.hist(d["speed"], bins = "auto")
plt.xlabel(xlabel)
plt.ylabel("number of occurrences")
plt.figure()
print("speed, when positive:")
- plot_distr_funct(speed[speed >= 0], label = "positive")
+ plot_distr_funct(d["speed"][d["speed"] >= 0], label = "positive")
print("- speed, when negative:")
- plot_distr_funct(- speed[speed < 0], label = "negative")
+ plot_distr_funct(- d["speed"][d["speed"] < 0], label = "negative")
plt.xlabel(xlabel)
plt.legend()
xlabel = "equivalent radius of outermost contour, in km"
plt.figure()
- plt.hist(rad_outer, bins = "auto")
+ plt.hist(d["rad_outer"], bins = "auto")
plt.xlabel(xlabel)
plt.ylabel("number of occurrences")
plt.figure()
print(xlabel + ":")
- plot_distr_funct(rad_outer)
+ plot_distr_funct(d["rad_outer"])
plt.xlabel(xlabel)
plt.ylabel("distribution function")
xlabel = "equivalent radius of max-speed contour, in km"
plt.figure()
- plt.hist(rad_speed, bins = "auto")
+ plt.hist(d["rad_speed"], bins = "auto")
plt.xlabel(xlabel)
plt.ylabel("number of occurrences")
plt.figure()
print(xlabel + ":")
- plot_distr_funct(rad_speed)
+ plot_distr_funct(d["rad_speed"])
plt.xlabel(xlabel)
plt.ylabel("distribution function")
xlabel = "amplitude of outermost contour, in m"
plt.figure()
- plt.hist(amp_outer, bins = "auto")
+ plt.hist(d["amp_outer"], bins = "auto")
plt.xlabel(xlabel)
plt.ylabel("number of occurrences")
plt.figure()
print("amplitude of outermost contour, anticyclones:")
- plot_distr_funct(amp_outer[amp_outer >= 0], label = "anticyclones")
+ plot_distr_funct(d["amp_outer"][d["amp_outer"] >= 0],
+ label = "anticyclones")
print("amplitude of outermost contour, cyclones:")
- plot_distr_funct(- amp_outer[amp_outer < 0], label = "cyclones")
+ plot_distr_funct(- d["amp_outer"][d["amp_outer"] < 0], label = "cyclones")
plt.xlabel(xlabel)
plt.ylabel("distribution function")
plt.legend()
@@ -121,15 +123,16 @@ if __name__ == "__main__":
xlabel = "amplitude of max-speed contour, in m"
plt.figure()
- plt.hist(amp_speed, bins = "auto")
+ plt.hist(d["amp_speed"], bins = "auto")
plt.xlabel(xlabel)
plt.ylabel("number of occurrences")
plt.figure()
print("amplitude of max-speed contour, anticyclones:")
- plot_distr_funct(amp_speed[amp_speed >= 0], label = "anticyclones")
+ plot_distr_funct(d["amp_speed"][d["amp_speed"] >= 0],
+ label = "anticyclones")
print("amplitude of max-speed contour, cyclones:")
- plot_distr_funct(- amp_speed[amp_speed < 0], label = "cyclones")
+ plot_distr_funct(- d["amp_speed"][d["amp_speed"] < 0], label = "cyclones")
plt.xlabel(xlabel)
plt.ylabel("distribution function")
plt.legend()
diff --git a/Analysis/plot_snapshot.py b/Analysis/plot_snapshot.py
index ec5ded35..8ef3ff0e 100755
--- a/Analysis/plot_snapshot.py
+++ b/Analysis/plot_snapshot.py
@@ -4,10 +4,11 @@
Red for anti-cyclones, blue for cyclones. Empty circles for extrema
with a valid outermost contour, empty squares for other
-extrema. Squares for a well-defined but invalid outermost
-contour. Crosses for a valid outermost contour but with no distinct
-max-speed contour. Filled circles for a valid outermost contour with a
-distinct max-speed contour.
+extrema. Squares on outermost contour for a well-defined but invalid
+outermost contour. Crosses on outermost contour for a valid outermost
+contour but with no distinct max-speed contour. Filled circles on
+outermost contour and max-speed contour for a valid outermost contour
+with a distinct max-speed contour.
This scripts takes about 30 s of CPU for a 90° by 90° window, about 4
mn for a 180° by 180° window.
@@ -39,15 +40,27 @@ parser.add_argument("-e", "--end_filename", help = "add character string after "
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.")
+ w = shapefile.Writer(shapeType = shapefile.NULL)
+ m_s_iterShapes = itertools.repeat(w)
+else:
+ m_s_iterShapes = reader_m_s.iterShapes()
+
with netCDF4.Dataset("h.nc") as f:
longitude = f.variables["lon"][:]
latitude = f.variables["lat"][:]
if args.window:
- llcrnrlon = float(input("llcrnrlon = ? "))
- llcrnrlat = float(input("llcrnrlat = ? "))
- urcrnrlon = float(input("urcrnrlon = ? "))
- urcrnrlat = float(input("urcrnrlat = ? "))
+ 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")
@@ -64,6 +77,7 @@ else:
plt.figure()
src_crs = ccrs.PlateCarree()
projection = ccrs.PlateCarree()
+##projection = ccrs.NorthPolarStereo()
ax = plt.axes(projection = projection)
if args.grid:
@@ -71,19 +85,6 @@ if args.grid:
ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
marker = "+", color = "gray", linestyle = "None")
-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.")
- w = shapefile.Writer(shapeType = shapefile.NULL)
- m_s_iterShapes = itertools.repeat(w)
-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]
@@ -103,7 +104,7 @@ for shape_rec_extr, shape_outer, shape_m_s \
lines = ax.plot(points[0], points[1], markersize = 10, color = color,
fillstyle = "none", transform = src_crs)
- if shape_rec_extr.record.valid != 0:
+ if shape_rec_extr.record.valid == 1:
if args.light:
lines[0].set_marker("+")
else:
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index 3917c53f..5d86f94a 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -20,6 +20,8 @@
\usepackage{algorithm}
+\newcommand{\ud}{\mathrm{d}}
+
\renewcommand{\algorithmicfor}{\textbf{pour}}
\renewcommand{\algorithmicend}{\textbf{fin}}
\renewcommand{\algorithmicdo}{\textbf{faire}}
@@ -1277,6 +1279,55 @@ weight à u et v, c'est-à -dire qu'il faudrait conserver u, v Ã
plusieurs dates en dehors de get\_snpashot. Cela remettrait en
question tout l'algorithme.
+\subsection{spherical\_polyline\_area}
+
+La procédure de calcul des contours donne une suite de sommets. Les
+courbes reliant ces sommets ne sont en fait pas imposés
+intrinsèquement par cette procédure. Par exemple, il n'est ni plus ni
+moins juste de supposer que ces courbes sont des portions de grands
+cercles ou de supposer qu'elles donnent des segments en projection
+plate carrée. Néanmoins, dans le reste du programme, pour décider si
+un contour entoure un extremum, ou pour calculer l'intersection de
+deux tourbillons, nous supposons que les contours sont des polygones
+en $(\lambda, \phi)$ (c'est-Ã -dire que les courbes reliant les sommets
+donnent des segments en projection plate carrée).
+
+Un domaine $S$ quelconque à la surface de la Terre a pour aire :
+\begin{equation*}
+ a^2 \iint_S \cos \phi\ \ud \phi\ \ud \lambda
+ = a^2 \iint_S \ud \lambda\ \ud(\sin \phi)
+\end{equation*}
+Nous avons une procédure qui calcule la surface $\iint \ud x\ \ud y$
+d'un polygone en $(x, y)$. Nous calculons la surface du polygone en
+$(\lambda, \sin \phi)$ passant par les sommets du
+contour. C'est-Ã -dire que nous supposons que les courbes reliant les
+somments donnent des segments en projection conforme cylindrique.
+
+Rémi fait un autre calcul : la surface du polygone en
+$(\lambda \cos \phi, \phi)$ passant par les sommets du contour. Le
+changement de variable :
+\begin{align*}
+ & x = \lambda \cos \phi \\
+ & y = \phi
+\end{align*}
+a pour jacobien $\cos \phi$ donc on a bien :
+\begin{equation*}
+ \iint \cos \phi\ \ud \phi\ \ud \lambda = \iint \ud x\ \ud y
+\end{equation*}
+Les deux surfaces obtenues, celle du polygone en
+$(\lambda, \sin \phi)$ et celle du polygone en
+$(\lambda \cos \phi, \phi)$ sont très proches si les sommets sont très
+proches. C'est le cas puisque les sommets sont séparés au maximum de
+\np{0.25}° en longitude et latitude. On peut aussi le voir en
+considérant la formule de calcul de surface du polygone :
+\begin{align*}
+ S & = \frac{1}{2} \sum_i x_i (y_{i + 1} - y_{i - 1}) \\
+ & =
+ \frac{1}{2} \sum_i \lambda_i \cos \phi_i (\phi_{i + 1} - \phi_{i - 1}) \\
+ & \approx \frac{1}{2} \sum_i \lambda_i (\sin \phi_{i + 1} - \sin \phi_{i - 1})
+\end{align*}
+si les valeurs successives des $\phi_i$ sont proches.
+
\subsection{successive\_overlap}
j, k, scalaires entiers, données. flow : donnée-résultat. Cherche les
diff --git a/spherical_polygon_area.f b/spherical_polygon_area.f
index d10b3e74..9fe5c848 100644
--- a/spherical_polygon_area.f
+++ b/spherical_polygon_area.f
@@ -6,8 +6,8 @@ contains
pure real function spherical_polygon_area(p)
- ! Assuming p is a polygon in longitude, latitude, in rad, compute
- ! the area of a polygon in longitude, sin(latitude) with the same
+ ! Assuming p is a polygon in longitude, latitude, compute the area
+ ! of a polygon in longitude, sin(latitude) with the same
! vertices. Result in m2.
! Libraries:
@@ -15,7 +15,7 @@ contains
use spherical_polyline_area_m, only: spherical_polyline_area
- type(polygon), intent(in):: p
+ type(polygon), intent(in):: p ! in rad
! Local:
diff --git a/spherical_polyline_area.f b/spherical_polyline_area.f
index b9d5449b..5c00b84b 100644
--- a/spherical_polyline_area.f
+++ b/spherical_polyline_area.f
@@ -6,15 +6,15 @@ contains
pure real function spherical_polyline_area(p)
- ! Assuming p is a polyline in longitude, latitude, in rad, compute
- ! the area (positive) of a polyline in longitude, sin(latitude)
- ! with the same vertices. Result in m2.
+ ! Assuming p is a polyline in longitude, latitude, compute the
+ ! area (positive) of a polyline in longitude, sin(latitude) with
+ ! the same vertices. Result in m2.
! Libraries:
use contour_531, only: polyline
use geometry, only: polygon_area_2d
- type(polyline), intent(in):: p ! should be closed
+ type(polyline), intent(in):: p ! Should be closed. In rad.
! Local:
--
GitLab