From c7cd342d65e37ef8a51e97397fa02eeed0fb17f6 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 8 Aug 2019 19:27:02 +0200
Subject: [PATCH] Add script `eddy_dump.py`.
Take into account new policy in `filter.py`: the names of the three
shapefiles are always the same and form a unit in a snaphot directory.
Rename `compil_prod_dir`to `build_dir` in JSON test files. Add a test
of `successive_overlap` with global grid, normal resolution.
---
Analysis/eddy_dump.py | 28 ++++++
Analysis/filter.py | 12 +--
Documentation_texfol/documentation.tex | 17 ++--
Tests/Input/successive_overlap_global_nml.txt | 2 +
Tests/long_tests.json | 8 +-
Tests/short_tests.json | 88 ++++++++++---------
extraction_eddies.f90 | 3 +
get_1_outerm.f90 | 6 +-
8 files changed, 105 insertions(+), 59 deletions(-)
create mode 100755 Analysis/eddy_dump.py
create mode 100644 Tests/Input/successive_overlap_global_nml.txt
diff --git a/Analysis/eddy_dump.py b/Analysis/eddy_dump.py
new file mode 100755
index 00000000..7e2eb26e
--- /dev/null
+++ b/Analysis/eddy_dump.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python3
+
+import argparse
+import shapefile
+from os import path
+import pprint
+
+parser = argparse.ArgumentParser()
+parser.add_argument("directory", help = "containing the three shapefiles")
+args = parser.parse_args()
+
+eddy_index = input("eddy_index = ? ")
+eddy_index = int(eddy_index)
+ishape = eddy_index - 1
+
+for basename in ["extremum", "outermost_contour", "max_speed_contour"]:
+ filename = path.join(args.directory, basename)
+
+ with shapefile.Reader(filename) as f:
+ shapeRec = f.shapeRecord(ishape)
+ print("\n", filename, ":")
+
+ if shapeRec.shape.shapeType == shapefile.NULL:
+ dct = shapeRec.record.as_dict()
+ pprint.pprint(dct)
+ print("null shape")
+ else:
+ pprint.pprint(shapeRec.__geo_interface__)
diff --git a/Analysis/filter.py b/Analysis/filter.py
index 4ed0faa0..0178d26b 100755
--- a/Analysis/filter.py
+++ b/Analysis/filter.py
@@ -2,12 +2,12 @@
import shapefile
-with shapefile.Reader("extremum_1") as extremum, \
- shapefile.Reader("outermost_contour_1") as outermost_cont, \
- shapefile.Reader("max_speed_contour_1") as max_speed_cont, \
- shapefile.Writer("extremum_1_filtered") as extremum_filt, \
- shapefile.Writer("outermost_contour_1_filtered") as outermost_cont_filt, \
- shapefile.Writer("max_speed_contour_1_filtered") as max_speed_cont_filt:
+with shapefile.Reader("extremum") as extremum, \
+ shapefile.Reader("outermost_contour") as outermost_cont, \
+ shapefile.Reader("max_speed_contour") as max_speed_cont, \
+ shapefile.Writer("extremum_filtered") as extremum_filt, \
+ shapefile.Writer("outermost_contour_filtered") as outermost_cont_filt, \
+ shapefile.Writer("max_speed_contour_filtered") as max_speed_cont_filt:
extremum_filt.fields = extremum.fields[1:]
outermost_cont_filt.fields = outermost_cont.fields[1:]
max_speed_cont_filt.fields = max_speed_cont.fields[1:]
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index df34d88e..8fb48e70 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -62,12 +62,17 @@ programme ne traite qu'une date.
\subsection{Cas pathologiques}
\verb+out_cont+ est vide si et seulement s'il n'y a
-pas de bon contour à \verb+innermost_level+. Le champ speed de
-extremum.dbf vaut \verb+missing_speed+ si et seulement si radius4 == 0
-ou \verb+speed_outerm+ dans \verb+set_max_speed+ est NaN. Si radius4
-$\ge$ 2 et speed dans extremum.dbf vaut \verb+missing_speed+ alors
-\verb+speed_outerm+ dans \verb+set_max_speed+ est NaN,
-\verb+e%speed_cont+ est
+pas de bon contour à \verb+innermost_level+, ce qui peut en
+particulier se produire si \verb+innermost_level+ correspond Ã
+l'amplitude minimale. Même si un contour extérieur est trouvé (avec
+une amplitude supérieure à l'amplitude minimale par construction), il
+peut ne pas être valide à cause du critère de surface minimale.
+
+Le champ speed de extremum.dbf vaut \verb+missing_speed+ si et
+seulement si radius4 == 0 ou \verb+speed_outerm+ dans
+\verb+set_max_speed+ est NaN. Si radius4 $\ge$ 2 et speed dans
+extremum.dbf vaut \verb+missing_speed+ alors \verb+speed_outerm+ dans
+\verb+set_max_speed+ est NaN, \verb+e%speed_cont+ est
non vide dans \verb+max_speed_contour.shp+ et \verb+mean_speed+ sur
\verb+e%speed_cont+ a donné un NaN. En d'autres termes, les deux
appels à \verb+mean_speed+ ont tous les deux produit un
diff --git a/Tests/Input/successive_overlap_global_nml.txt b/Tests/Input/successive_overlap_global_nml.txt
new file mode 100644
index 00000000..bcbe3d6d
--- /dev/null
+++ b/Tests/Input/successive_overlap_global_nml.txt
@@ -0,0 +1,2 @@
+&MAIN_NML CORNER_DEG=0.125,-89.875, NLON =1440 NLAT =720/
+
diff --git a/Tests/long_tests.json b/Tests/long_tests.json
index 8365bf26..4bd4aede 100644
--- a/Tests/long_tests.json
+++ b/Tests/long_tests.json
@@ -1,12 +1,12 @@
[
{
- "args": ["$compil_prod_dir/test_set_all_outerm",
+ "args": ["$build_dir/test_set_all_outerm",
"$input_dir/h_region_3.nc"],
"title": "Set_all_outerm",
"input": "&main_nml min_amp = 0.001/\n"
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_3",
"required": [["$input_dir/h_region_3.nc", "h.nc"],
["$input_dir/uv_region_3.nc", "uv.nc"]],
@@ -14,7 +14,7 @@
"description": "Larger region, 120 x 120. Includes degenerate extrema."
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_3_min",
"required": [["$input_dir/h_region_3.nc", "h.nc"],
["$input_dir/uv_region_3.nc", "uv.nc"]],
@@ -22,7 +22,7 @@
"description": "Same as Extraction_eddies_region_3 except with 1 mm minimum amplitude."
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_5",
"required": [["$input_dir/h_region_5.nc", "h.nc"],
["$input_dir/uv_region_5.nc", "uv.nc"]],
diff --git a/Tests/short_tests.json b/Tests/short_tests.json
index 9c0ff18a..5f100604 100644
--- a/Tests/short_tests.json
+++ b/Tests/short_tests.json
@@ -1,7 +1,7 @@
[
{
"stdin_filename" : "$input_dir/good_contour.txt",
- "args" : ["$compil_prod_dir/test_good_contour",
+ "args" : ["$build_dir/test_good_contour",
"$input_dir/example.nc"],
"title" : "Good_contour",
"required": [["$input_dir/outside_points_1.csv", "outside_points.csv"]],
@@ -10,7 +10,7 @@
},
{
"stdin_filename" : "$input_dir/good_contour_2.txt",
- "args" : ["$compil_prod_dir/test_good_contour",
+ "args" : ["$build_dir/test_good_contour",
"$input_dir/example.nc"],
"title" : "Good_contour_2",
"required": [["$input_dir/outside_points_2.csv", "outside_points.csv"]],
@@ -19,26 +19,26 @@
},
{
"stdin_filename" : "$input_dir/no_good_contour.txt",
- "args" : ["$compil_prod_dir/test_good_contour",
+ "args" : ["$build_dir/test_good_contour",
"$input_dir/example.nc"],
"title" : "No_good_contour",
"required": [["$input_dir/outside_points_1.csv", "outside_points.csv"]]
},
{
- "args" : ["$compil_prod_dir/test_local_extrema",
+ "args" : ["$build_dir/test_local_extrema",
"$input_dir/h_region_1.nc"],
"title" : "Local_extrema",
"input" : "f"
},
{
- "args" : ["$compil_prod_dir/test_local_extrema",
+ "args" : ["$build_dir/test_local_extrema",
"$input_dir/h_region_2.nc"],
"title" : "Local_extrema_larger",
"input" : "f",
"description": "With a larger region than in test Local_extrema."
},
{
- "args" : "$compil_prod_dir/test_get_1_outerm",
+ "args" : "$build_dir/test_get_1_outerm",
"required":
[
["$input_dir/h_region_1.nc", "h.nc"],
@@ -48,7 +48,7 @@
"input" : "&main_nml /\n"
},
{
- "args" : "$compil_prod_dir/test_get_1_outerm",
+ "args" : "$build_dir/test_get_1_outerm",
"required":
[
["$input_dir/h_region_1.nc", "h.nc"],
@@ -59,7 +59,7 @@
"description": "Assume insufficient amplitude for extrema 2 and 8. Even if extremum 8 does not have a particularly small amplitude, we just want to see the target extremum, 6, grow its outermost contour compared to the case without minimum amplitude."
},
{
- "args" : "$compil_prod_dir/test_get_1_outerm",
+ "args" : "$build_dir/test_get_1_outerm",
"required": [
["$input_dir/h_region_1.nc", "h.nc"],
[
@@ -72,20 +72,20 @@
"description": "Negative value for extremum 2."
},
{
- "args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
+ "args" : ["$build_dir/test_max_speed_contour_ssh",
"$input_dir/h_region_1.nc", "$input_dir/uv_region_1.nc"],
"title" : "Max_speed_contour_ssh",
"input" : "&main_nml /\n"
},
{
- "args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
+ "args" : ["$build_dir/test_max_speed_contour_ssh",
"$input_dir/h_region_1.nc", "$input_dir/uv_region_1.nc"],
"title" : "Max_speed_contour_ssh_north",
"input" : "&main_nml IND_EXTR= 4,14/\n",
"description": "direction = 2"
},
{
- "args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
+ "args" : ["$build_dir/test_max_speed_contour_ssh",
"$input_dir/huv_2015_11_29.nc",
"$input_dir/huv_2015_11_29.nc"],
"title" : "Max_speed_contour_ssh_missing",
@@ -93,7 +93,7 @@
"description": "Missing value in speed."
},
{
- "args" : "$compil_prod_dir/test_mean_speed",
+ "args" : "$build_dir/test_mean_speed",
"required": [["$input_dir/uv_region_1.nc", "uv.nc"],
["$tests_old_dir/Get_1_outerm/test_get_1_outerm.shp",
"contour.shp"],
@@ -105,7 +105,7 @@
"input" : "&main_nml /\n"
},
{
- "args" : "$compil_prod_dir/test_mean_speed",
+ "args" : "$build_dir/test_mean_speed",
"required": [["$input_dir/uv_region_1.nc","uv.nc"],
["$input_dir/outermost_contour_alt.shp", "contour.shp"],
["$input_dir/outermost_contour_alt.dbf", "contour.dbf"],
@@ -114,7 +114,7 @@
"input" : "&main_nml /\n"
},
{
- "args": "$compil_prod_dir/test_set_max_speed",
+ "args": "$build_dir/test_set_max_speed",
"title": "Set_max_speed",
"required": [["$input_dir/h_outermost.nc", "h.nc"],
["$input_dir/uv_outermost.nc", "uv.nc"],
@@ -128,7 +128,7 @@
"input" : "&main_nml /\n"
},
{
- "args": "$compil_prod_dir/test_set_max_speed",
+ "args": "$build_dir/test_set_max_speed",
"title": "Set_max_speed_noise",
"required":
[
@@ -140,14 +140,14 @@
"input": "&MAIN_NML IND_TARG_EXTR=19,11 /\n"
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_1",
"required": [["$input_dir/h_region_1.nc", "h.nc"],
["$input_dir/uv_region_1.nc", "uv.nc"]],
"input": "&main_nml min_amp = 0./\n"
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_1_noise",
"required": [["$input_dir/h_region_1.nc", "h.nc"],
["$input_dir/uv_region_1.nc", "uv.nc"]],
@@ -155,33 +155,33 @@
"description" : "Same as Extraction_eddies_region_1 but with non-zero minimal amplitude."
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_2",
"required": [["$input_dir/h_region_2.nc", "h.nc"],
["$input_dir/uv_region_2.nc", "uv.nc"]],
"input": "&main_nml min_amp = 0./\n"
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_2_noise",
"required": [["$input_dir/h_region_2.nc", "h.nc"],
["$input_dir/uv_region_2.nc", "uv.nc"]],
"input": "&main_nml /\n"
},
{
- "args": ["$compil_prod_dir/test_inside_4",
+ "args": ["$build_dir/test_inside_4",
"$input_dir/outermost_eddy_5"],
"title": "Inside_4_true",
"stdin_filename": "$input_dir/inside_4_true_nml.txt"
},
{
- "args": ["$compil_prod_dir/test_inside_4",
+ "args": ["$build_dir/test_inside_4",
"$input_dir/outermost_eddy_5"],
"title": "Inside_4_false",
"stdin_filename": "$input_dir/inside_4_false_nml.txt"
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_region_4",
"required":
[
@@ -191,44 +191,44 @@
"input": "&main_nml /\n",
"description": "Part of the domain has missing values."
},
- {"args": "$compil_prod_dir/test_weight", "title": "Weight"},
+ {"args": "$build_dir/test_weight", "title": "Weight"},
{
- "args": "$compil_prod_dir/test_spher_polyline_area",
+ "args": "$build_dir/test_spher_polyline_area",
"title": "Spher_polyline_area"
},
{
- "args": ["$compil_prod_dir/test_spher_polygon_area",
+ "args": ["$build_dir/test_spher_polygon_area",
"$input_dir/triangle"],
"title": "Spher_polygon_area"
},
{
- "args": ["$compil_prod_dir/test_spher_polygon_area",
+ "args": ["$build_dir/test_spher_polygon_area",
"$input_dir/tri_square_hole"],
"title": "Area_multi_polygon",
"description": "Area of a multipolygon, with a hole in one polygon."
},
{
- "args": "$compil_prod_dir/test_read_eddy.sh",
+ "args": "$build_dir/test_read_eddy.sh",
"required":
[["$input_dir/Region_4_2006_01_01/Snapshot", "Snapshot_old"]],
"title": "Read_eddy"
},
{
- "args": "$compil_prod_dir/test_read_snapshot.sh",
+ "args": "$build_dir/test_read_snapshot.sh",
"required":
[["$input_dir/Region_4_2006_01_01/Snapshot", "Snapshot_old"]],
"title": "Read_snapshot",
"stdin_filename": "$input_dir/read_snapshot_nml.txt"
},
{
- "args" : ["$compil_prod_dir/test_local_extrema",
+ "args" : ["$build_dir/test_local_extrema",
"$input_dir/Region_4_2006_01_01/h.nc"],
"title" : "Local_extrema_missing",
"description": "With input file containing missing values.",
"input" : "f"
},
{
- "args" : ["$compil_prod_dir/test_successive_overlap",
+ "args" : ["$build_dir/test_successive_overlap",
"$input_dir/Region_4_2006_01_01/Snapshot",
"$input_dir/Region_4_2006_01_01/Snapshot"],
"title" : "Successive_overlap",
@@ -236,7 +236,7 @@
"stdin_filename": "$input_dir/successive_overlap_nml.txt"
},
{
- "args" : "$compil_prod_dir/test_nearby_extr",
+ "args" : "$build_dir/test_nearby_extr",
"title" : "Nearby_extr",
"required":
[
@@ -245,7 +245,7 @@
]
},
{
- "args" : ["$compil_prod_dir/test_successive_overlap",
+ "args" : ["$build_dir/test_successive_overlap",
"$tests_old_dir/Extraction_eddies_region_5/Snapshot",
"$tests_old_dir/Extraction_eddies_region_5/Snapshot"],
"title" : "Successive_overlap_region_5_one_date",
@@ -253,28 +253,28 @@
"stdin_filename": "$input_dir/successive_overlap_region_5_nml.txt"
},
{
- "args" : ["$compil_prod_dir/test_local_extrema",
+ "args" : ["$build_dir/test_local_extrema",
"$input_dir/h_region_1.nc"],
"title" : "Local_extrema_periodic",
"input" : "t",
"description": "Same as Local_extrema but with periodicity. The data is actually regional so there is a discontinuity in the field and the longitude grid assuming periodicity is not regular."
},
{
- "args" : ["$compil_prod_dir/test_local_extrema",
+ "args" : ["$build_dir/test_local_extrema",
"$input_dir/h_2006_01_01_coarse.nc"],
"title" : "Local_extrema_periodic_2",
"input" : "t",
"description": "The data is really global so there is no discontinuity in the field and the longitude grid assuming periodicity is regular."
},
{
- "args": ["$compil_prod_dir/test_set_all_outerm",
+ "args": ["$build_dir/test_set_all_outerm",
"$input_dir/h_2006_01_01_coarse.nc"],
"title": "Set_all_outerm_periodic",
"input": "&main_nml MAX_RADIUS = 4, 4/\n",
"description": "test_set_all_outerm with periodicity."
},
{
- "args": "$compil_prod_dir/extraction_eddies.sh",
+ "args": "$build_dir/extraction_eddies.sh",
"title": "Extraction_eddies_periodic",
"required": [["$input_dir/h_2006_01_01_coarse.nc", "h.nc"],
["$input_dir/uv_2006_01_01_coarse.nc", "uv.nc"]],
@@ -282,7 +282,7 @@
"description": "Periodic domain."
},
{
- "args" : ["$compil_prod_dir/test_successive_overlap",
+ "args" : ["$build_dir/test_successive_overlap",
"$tests_old_dir/Extraction_eddies_periodic/Snapshot",
"$tests_old_dir/Extraction_eddies_periodic/Snapshot"],
"title" : "Successive_overlap_periodic",
@@ -290,7 +290,7 @@
"stdin_filename": "$input_dir/successive_overlap_periodic_nml.txt"
},
{
- "args" : ["$compil_prod_dir/test_successive_overlap",
+ "args" : ["$build_dir/test_successive_overlap",
"$input_dir/Region_4_2006_01_01/Snapshot",
"$input_dir/Region_4_2006_01_02/Snapshot"],
"title" : "Successive_overlap_different_snapshots",
@@ -298,7 +298,7 @@
"stdin_filename": "$input_dir/successive_overlap_nml.txt"
},
{
- "args" : ["$compil_prod_dir/test_successive_overlap",
+ "args" : ["$build_dir/test_successive_overlap",
"$tests_old_dir/Extraction_eddies_region_2_noise/Snapshot",
"$input_dir/Region_2_2006_01_02/Snapshot"],
"title" : "Successive_overlap_region_2",
@@ -306,11 +306,19 @@
"stdin_filename": "$input_dir/successive_overlap_region_2_nml.txt"
},
{
- "args" : ["$compil_prod_dir/test_successive_overlap",
+ "args" : ["$build_dir/test_successive_overlap",
"$tests_old_dir/Extraction_eddies_region_5/Snapshot",
"$input_dir/Snapshot_region_5_2006_01_02"],
"title" : "Successive_overlap_region_5",
"description": "Same as Successive_overlap_region_2, but with a larger region. We get some merging and splitting.",
"stdin_filename": "$input_dir/successive_overlap_region_5_nml.txt"
+ },
+ {
+ "args" : ["$build_dir/test_successive_overlap",
+ "$large_input_dir/Snapshot_2006_01_01",
+ "$large_input_dir/Snapshot_2006_01_02"],
+ "title" : "Successive_overlap_global",
+ "description": "Global grid, normal 0.25° resolution.",
+ "stdin_filename": "$input_dir/successive_overlap_global_nml.txt"
}
]
diff --git a/extraction_eddies.f90 b/extraction_eddies.f90
index 5f04e0a6..e7ef1363 100644
--- a/extraction_eddies.f90
+++ b/extraction_eddies.f90
@@ -77,6 +77,7 @@ program extraction_eddies
call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
"extraction_eddies: not enough IEEE support")
+ ! We need this before looking at the input files:
write(unit = error_unit, nml = main_nml)
write(unit = error_unit, fmt = *) "Enter namelist main_nml."
read(unit = *, nml = main_nml)
@@ -88,8 +89,10 @@ program extraction_eddies
call nf95_get_att(ncid, varid, name = "units", values = time_unit)
call assert(time_unit == "days since 1950-01-01 00:00:00" &
.or. time_unit == "days since 1950-1-1 00:00:00", "main: bad time unit")
+
! We are assuming there is a single date in the input file:
call nf95_get_var(ncid, varid, time)
+
k = nint(time)
call assert(abs(time - k) < 0.1, "main: bad time value")
diff --git a/get_1_outerm.f90 b/get_1_outerm.f90
index 5017fbf1..01657288 100644
--- a/get_1_outerm.f90
+++ b/get_1_outerm.f90
@@ -82,7 +82,7 @@ contains
get_1_outerm%polyline = tentative_contour
get_1_outerm%ssh = level_try
else
- ! {no good contour at maxval or minval}
+ ! Assertion: {no good contour at maxval or minval}
level_bad = level_try
do while (abs(level_bad - level_good) > precision)
@@ -96,9 +96,9 @@ contains
level_bad = level_try
end if
- ! {get_1_outerm%polyline == good_contour(. . .,
+ ! Assertion: get_1_outerm%polyline == good_contour(. . .,
! level_good) and get_1_outerm%%n_points /= 0 and
- ! no good contour at level_bad}
+ ! no good contour at level_bad
end do
get_1_outerm%ssh = level_good
--
GitLab