From a6ab40be580f6ad0d388b1b164c92607d5e502e0 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Tue, 4 Dec 2018 15:47:48 +0100
Subject: [PATCH] For clarity, in procedure set_max_speed, do not modify
speed_outerm after its initial computation and show more directly the
consequence of a NaN value. The drawback is the repetition of one line of
code. As for performance, one less test is performed in the case of a NaN
value.
In script stat.py, add stats for cases when speed is missing.
---
Documentation_texfol/documentation.tex | 9 ++++++++
Tests/stat.py | 25 ++++++++++++++++++++-
get_snapshot.f | 2 ++
set_all_outerm.f | 3 ++-
set_max_speed.f | 30 +++++++++++++++-----------
5 files changed, 55 insertions(+), 14 deletions(-)
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index de818de5..ad064dbf 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -111,6 +111,15 @@ Cf. figure (\ref{fig:call_graph}).
\label{fig:call_graph}
\end{figure}
+Cas pathologiques. out\_cont est vide si et seulement s'il n'y a pas
+de bon contour à innermost\_level. Le champ speed de extremum.dbf vaut
+missing\_speed si et seulement si radius4 == 0 ou speed\_outerm dans
+set\_max\_speed est NaN. Si radius4 $\ge$ 2 et speed dans extremum.dbf
+vaut missing\_speed alors speed\_outerm dans set\_max\_speed est NaN,
+e\%speed\_cont est non vide dans max\_speed\_contour.shp et
+mean\_speed sur e\%speed\_cont a donné un NaN. En d'autres termes, les
+deux appels à mean\_speed ont tous les deux produit un NaN.
+
Cas où les composantes out\_cont ou speed\_cont d'un tourbillon sont
vides, c'est-à-dire égales à null\_ssh\_contour(). On peut distinguer
les cas du tableau (\ref{tab:null_ssh_contour}).
diff --git a/Tests/stat.py b/Tests/stat.py
index f12b46d7..fb2de328 100755
--- a/Tests/stat.py
+++ b/Tests/stat.py
@@ -10,6 +10,7 @@ n_radius4 = np.zeros(3, dtype = int)
n_outer = 0
n_points_valid = 0
n_points = 0
+n_missing_speed = np.zeros(3, dtype = int)
with shapefile.Reader("extremum_1") as extremum, \
shapefile.Reader("outermost_contour_1") as outermost_contour:
@@ -35,23 +36,45 @@ with shapefile.Reader("extremum_1") as extremum, \
else:
n_radius4[2] += 1
-print("number of extrema of SSH:", len(extremum))
+ if rec_extr.speed == 1e4:
+ if shape_rec_outer.record.radius4 == 0:
+ n_missing_speed[0] += 1
+ elif shape_rec_outer.record.radius4 == 1:
+ n_missing_speed[1] += 1
+ else:
+ n_missing_speed[2] += 1
+
+print("number of extrema of SSH:", n_extr)
print("number of minima of SSH:", n_minima)
+print("number of maxima of SSH:", n_extr - n_minima)
+print()
+
print("number of valid eddies (outermost contour found with sufficient area):",
n_valid)
print("number of valid cyclones:", n_valid_cycl)
+print("number of valid anticyclones:", n_valid - n_valid_cycl)
+print()
+
print("number of extrema with radius4 == 0 (no outermost contour or "
"insufficient area):", n_radius4[0])
print("number of extrema with radius4 == 1 (valid eddies with radius4 == 1):",
n_radius4[1])
print("number of extrema with radius4 >= 2 (valid eddies with radius4 >= 2):",
n_radius4[2])
+print()
+
print("number of outermost contours (outermost contour found, but may be with "
"sufficient area or not):", n_outer)
print("average number of points of outermost contour per valid eddy:",
n_points_valid / n_valid)
print("average number of points of outermost contour per extremum (valid eddy "
"or not):", n_points / n_extr)
+print()
+
+print("number of missing speed values with radius4 == 0:", n_missing_speed[0])
+print("number of missing speed values with radius4 == 1:", n_missing_speed[1])
+print("number of missing speed values with radius4 >= 2:", n_missing_speed[2])
+print()
n_max_speed = 0
n_points_valid = 0
diff --git a/get_snapshot.f b/get_snapshot.f
index f2bb0a5b..fa8261b7 100644
--- a/get_snapshot.f
+++ b/get_snapshot.f
@@ -113,11 +113,13 @@ contains
outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
llc(2):urc(2)), s%list_vis, i)
+
if (periodic) outside_points(1, :) = outside_points(1, :) &
+ ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
* twopi
! (Shift the longitude of each point to a value close to
! the longitude of the target extremum.)
+
call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
outside_points, ssh(llc(1):urc(1), llc(2):urc(2)), &
u(llc(1):urc(1), llc(2):urc(2)), &
diff --git a/set_all_outerm.f b/set_all_outerm.f
index 194728d6..178c336f 100644
--- a/set_all_outerm.f
+++ b/set_all_outerm.f
@@ -22,7 +22,8 @@ contains
type(snapshot), intent(out):: s
! Specifically: define everything in s except
- ! s%list_vis%speed_cont, s%list_vis%max_speed and s%number_eddies.
+ ! s%list_vis%speed_cont, s%list_vis%max_speed, s%list_vis%radius4
+ ! and s%number_eddies.
real, intent(in):: min_amp ! minimum amplitude of ssh, between
! extremum and outermost contour, in m
diff --git a/set_max_speed.f b/set_max_speed.f
index 47b435eb..5ba8879d 100644
--- a/set_max_speed.f
+++ b/set_max_speed.f
@@ -8,9 +8,10 @@ contains
corner, step)
! Defines the components speed_cont, max_speed and radius4 of
- ! argument e. On return, speed_cont may be a null ssh contour but
- ! max_speed is never set to missing_speed. The value of radius4
- ! set here is always >= 1.
+ ! argument e. On return, e%speed_cont may be a null ssh contour
+ ! but e%max_speed is never set to missing_speed. However,
+ ! e%max_speed might be NaN if speed_outerm is NaN (and any value
+ ! of e%radius4). The value of e%radius4 set here is always >= 1.
! The length of the interval of longitudes of the domain, step(1)
! * (size(ssh, 1) - 1), should be < 180 degrees, so that the
@@ -24,7 +25,7 @@ contains
! Lbraries:
use contour_531, only: convert_to_ind, null_polyline
use geometry, only: polygon_point_dist_2d
-
+
use derived_types, only: eddy, null_ssh_contour
use good_contour_m, only: good_contour
use inside_4_m, only: inside_4
@@ -85,16 +86,21 @@ contains
e%max_speed = mean_speed(u, v, e%speed_cont%polyline, e%coord_extr, &
corner, step)
- if (IEEE_IS_NAN(speed_outerm)) speed_outerm = 0.
- ! (This may happen when the outermost contour is near land.)
-
- if (abs(speed_outerm) > abs(e%max_speed)) then
- ! Abandon the contour coming from max_speed_contour_ssh:
- e%speed_cont = null_ssh_contour()
- e%max_speed = speed_outerm
- else
+ if (IEEE_IS_NAN(speed_outerm)) then
+ ! This may happen when the outermost contour is near land.
! Stick to the contour coming from max_speed_contour_ssh.
e%speed_cont%area = spherical_polyline_area(e%speed_cont%polyline)
+ else
+ ! Note the following test should raise an invalid exception
+ ! if e%max_speed is NaN.
+ if (abs(speed_outerm) > abs(e%max_speed)) then
+ ! Abandon the contour coming from max_speed_contour_ssh:
+ e%speed_cont = null_ssh_contour()
+ e%max_speed = speed_outerm
+ else
+ ! Stick to the contour coming from max_speed_contour_ssh.
+ e%speed_cont%area = spherical_polyline_area(e%speed_cont%polyline)
+ end if
end if
end if
--
GitLab