From a33ecf3eda9dc116e0b0cd802a27ce22ebb43363 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 17 May 2018 16:32:24 +0200
Subject: [PATCH] Add procedure successive_overlap and corresponding test. The
libraries GPC_F and GPC must be used.
---
GNUmakefile | 10 +++-
Tests/Stdin/successive_overlap_nml.txt | 5 ++
Tests/read_snapshot.f | 10 +++-
Tests/test_successive_overlap.f | 57 ++++++++++++++++++
Tests/tests.json | 36 ++++++++++++
depend.mk | 2 +
derived_types.f | 12 ++--
successive_overlap.f | 81 ++++++++++++++++++++++++++
8 files changed, 202 insertions(+), 11 deletions(-)
create mode 100644 Tests/Stdin/successive_overlap_nml.txt
create mode 100644 Tests/test_successive_overlap.f
create mode 100644 successive_overlap.f
diff --git a/GNUmakefile b/GNUmakefile
index 2c198fef..5aae2083 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -23,9 +23,11 @@ src_test_read_eddy = test_read_eddy.f derived_types.f init_shapefiles.f read_edd
src_test_read_snapshot = test_read_snapshot.f derived_types.f dispatch_snapshot.f init_shapefiles.f read_snapshot.f write_eddy.f send_snapshot.f read_eddy.f
-sources := $(sort ${src_test_local_extrema} ${src_test_get_1_outerm} ${src_test_set_max_speed} ${src_test_get_snapshot} ${src_test_set_all_outerm} ${src_test_weight} ${src_test_spherical_polygon_area} ${src_test_read_eddy} ${src_test_read_snapshot}) test_good_contour.f test_inside_4.f test_max_speed_contour_ssh.f test_mean_speed.f test_spherical_polyline_area.f
+src_test_successive_overlap = test_successive_overlap.f derived_types.f successive_overlap.f read_snapshot.f spherical_polygon_area.f spherical_polyline_area.f weight.f read_eddy.f
-lib_list = GPC_F contour_531 numer_rec_95 shapelib_03 netcdf95 geometry jumble netcdff fortrangis shp fortranc nr_util
+sources := $(sort ${src_test_local_extrema} ${src_test_get_1_outerm} ${src_test_set_max_speed} ${src_test_get_snapshot} ${src_test_set_all_outerm} ${src_test_weight} ${src_test_spherical_polygon_area} ${src_test_read_eddy} ${src_test_read_snapshot} ${src_test_successive_overlap}) test_good_contour.f test_inside_4.f test_max_speed_contour_ssh.f test_mean_speed.f test_spherical_polyline_area.f
+
+lib_list = GPC_F contour_531 numer_rec_95 GPC shapelib_03 netcdf95 geometry jumble netcdff fortrangis shp fortranc nr_util
# 2. Objects and executable files
@@ -38,10 +40,11 @@ obj_test_weight := $(src_test_weight:.f=.o)
obj_test_spherical_polygon_area := $(src_test_spherical_polygon_area:.f=.o)
obj_test_read_eddy := $(src_test_read_eddy:.f=.o)
obj_test_read_snapshot := $(src_test_read_snapshot:.f=.o)
+obj_test_successive_overlap := $(src_test_successive_overlap:.f=.o)
objects := $(sources:.f=.o)
-execut = test_good_contour test_inside_4 test_get_1_outerm test_local_extrema test_max_speed_contour_ssh test_mean_speed test_set_max_speed test_get_snapshot test_set_all_outerm test_weight test_spherical_polyline_area test_spherical_polygon_area test_read_eddy test_read_snapshot
+execut = test_good_contour test_inside_4 test_get_1_outerm test_local_extrema test_max_speed_contour_ssh test_mean_speed test_set_max_speed test_get_snapshot test_set_all_outerm test_weight test_spherical_polyline_area test_spherical_polygon_area test_read_eddy test_read_snapshot test_successive_overlap
# 3. Compiler-dependent part
@@ -80,6 +83,7 @@ test_spherical_polyline_area: spherical_polyline_area.o
test_spherical_polygon_area: ${obj_test_spherical_polygon_area}
test_read_eddy: ${obj_test_read_eddy}
test_read_snapshot: ${obj_test_read_snapshot}
+test_successive_overlap: ${obj_test_successive_overlap}
depend ${makefile_dir}/depend.mk:
makedepf90 -free -Wmissing -Wconfused $(addprefix -I, ${VPATH}) -nosrc $(addprefix -u , ${lib_list} shapelib netcdf intrinsic) ${sources} >${makefile_dir}/depend.mk
diff --git a/Tests/Stdin/successive_overlap_nml.txt b/Tests/Stdin/successive_overlap_nml.txt
new file mode 100644
index 00000000..f4f7c4fe
--- /dev/null
+++ b/Tests/Stdin/successive_overlap_nml.txt
@@ -0,0 +1,5 @@
+ &MAIN_NML
+ CORNER = 16.12500 , -38.87500 ,
+ NLON = 20,
+ NLAT = 20
+ /
diff --git a/Tests/read_snapshot.f b/Tests/read_snapshot.f
index dfee0586..3ac2cd43 100644
--- a/Tests/read_snapshot.f
+++ b/Tests/read_snapshot.f
@@ -11,7 +11,7 @@ contains
use contour_531, only: convert_to_ind
use derived_types, only: snapshot, eddy
- use nr_util, only: pi, assert, deg_to_rad
+ use nr_util, only: assert, deg_to_rad
use read_eddy_m, only: read_eddy
use shapelib, only: shpfileobject, shpgetinfo
@@ -30,7 +30,13 @@ contains
integer, intent(in):: nlon, nlat
! Local:
- integer ishape, k, i
+ integer ishape
+
+ integer k
+ ! Date index, not used. We could check that it has the same value
+ ! for all input shapes.
+
+ integer i ! eddy index
integer shapetype, dbffieldcount, dbfrecordcount
real(c_double) minbound(4), maxbound(4)
type(eddy) e
diff --git a/Tests/test_successive_overlap.f b/Tests/test_successive_overlap.f
new file mode 100644
index 00000000..e30e359f
--- /dev/null
+++ b/Tests/test_successive_overlap.f
@@ -0,0 +1,57 @@
+program test_successive_overlap
+
+ use, intrinsic:: ISO_FORTRAN_ENV
+
+ use derived_types, only: snapshot
+ use jumble, only: new_unit
+ use nr_util, only: deg_to_rad
+ use read_snapshot_m, only: read_snapshot
+ use shapelib, only: shpfileobject, shpclose
+ use shapelib_03, only: shp_open_03
+ use successive_overlap_m, only: successive_overlap
+
+ implicit none
+
+ integer unit_graph, i
+ type(snapshot) flow(2)
+ TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
+ TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
+ TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
+
+ real:: corner(2) = [0.125, - 59.875]
+ ! longitude and latitude of the corner of the whole grid, in degrees
+
+ integer:: nlon = 120, nlat = 120
+ namelist /main_nml/ corner, nlon, nlat
+
+ !-------------------------------------------------------------------------
+
+ write(unit = error_unit, nml = main_nml)
+ write(unit = error_unit, fmt = *) "Enter namelist main_nml."
+ read(unit = *, nml = main_nml)
+ write(unit = *, nml = main_nml)
+
+ call shp_open_03("extremum_1", pszaccess = "rb", hshp = hshp_extremum)
+ call shp_open_03("outermost_contour_1", pszaccess = "rb", &
+ hshp = hshp_outermost)
+ call shp_open_03("max_speed_contour_1", pszaccess = "rb", &
+ hshp = hshp_max_speed)
+ call read_snapshot(flow(1), hshp_extremum, hshp_outermost, hshp_max_speed, &
+ corner * deg_to_rad, nlon, nlat)
+ CALL shpclose(hshp_extremum)
+ CALL shpclose(hshp_outermost)
+ CALL shpclose(hshp_max_speed)
+
+ flow(2) = flow(1)
+ call new_unit(unit_graph)
+ open(unit_graph, file = "graph_1.txt", status = "replace", action = "write")
+ call successive_overlap(unit_graph, nlon, nlat, flow, j = 2, k = 2)
+ close(unit_graph)
+ print *, 'Created file "graph_1.txt".'
+
+ do i = 1, 2
+ print *, "flow(", i, ")%list_vis%delta_in:", flow(i)%list_vis%delta_in
+ print *, "flow(", i, ")%list_vis%delta_out:", flow(i)%list_vis%delta_out
+ end do
+
+end program test_successive_overlap
diff --git a/Tests/tests.json b/Tests/tests.json
index 792bf929..3cbc520e 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -219,5 +219,41 @@
"$input_dir/2006_01/h_region_4.nc"],
"title" : "Local_extrema_missing",
"description": "test_local_extrema with input file containing missing values."
+ },
+ {
+ "args" : "$compil_prod_dir/test_successive_overlap",
+ "title" : "Successive_overlap",
+ "description": "Overlap of a snapshot with itself.",
+ "required":
+ [
+ ["$input_dir/Test_read_eddy/extremum_1_old.dbf", "extremum_1.dbf"],
+ ["$input_dir/Test_read_eddy/extremum_1_old.shp", "extremum_1.shp"],
+ ["$input_dir/Test_read_eddy/extremum_1_old.shx", "extremum_1.shx"],
+ [
+ "$input_dir/Test_read_eddy/max_speed_contour_1_old.dbf",
+ "max_speed_contour_1.dbf"
+ ],
+ [
+ "$input_dir/Test_read_eddy/max_speed_contour_1_old.shp",
+ "max_speed_contour_1.shp"
+ ],
+ [
+ "$input_dir/Test_read_eddy/max_speed_contour_1_old.shx",
+ "max_speed_contour_1.shx"
+ ],
+ [
+ "$input_dir/Test_read_eddy/outermost_contour_1_old.dbf",
+ "outermost_contour_1.dbf"
+ ],
+ [
+ "$input_dir/Test_read_eddy/outermost_contour_1_old.shp",
+ "outermost_contour_1.shp"
+ ],
+ [
+ "$input_dir/Test_read_eddy/outermost_contour_1_old.shx",
+ "outermost_contour_1.shx"
+ ]
+ ],
+ "stdin": "$stdin_dir/successive_overlap_nml.txt"
}
]
diff --git a/depend.mk b/depend.mk
index 93ccfd2b..5fefc99a 100644
--- a/depend.mk
+++ b/depend.mk
@@ -8,6 +8,7 @@ send_snapshot.o : derived_types.o
set_all_outerm.o : local_extrema.o get_1_outerm.o derived_types.o
set_max_speed.o : spherical_polyline_area.o mean_speed.o max_speed_contour_ssh.o inside_4.o good_contour.o derived_types.o
spherical_polygon_area.o : spherical_polyline_area.o
+successive_overlap.o : weight.o spherical_polyline_area.o spherical_polygon_area.o derived_types.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
test_get_snapshot.o : init_shapefiles.o get_snapshot.o dispatch_snapshot.o derived_types.o
test_local_extrema.o : local_extrema.o
@@ -16,6 +17,7 @@ test_read_snapshot.o : read_snapshot.o init_shapefiles.o dispatch_snapshot.o der
test_set_all_outerm.o : set_all_outerm.o derived_types.o
test_set_max_speed.o : set_max_speed.o derived_types.o
test_spherical_polygon_area.o : spherical_polygon_area.o
+test_successive_overlap.o : successive_overlap.o read_snapshot.o derived_types.o
test_weight.o : weight.o derived_types.o
weight.o : derived_types.o
write_eddy.o : init_shapefiles.o derived_types.o
diff --git a/derived_types.f b/derived_types.f
index e5508d12..ae85197c 100644
--- a/derived_types.f
+++ b/derived_types.f
@@ -45,12 +45,12 @@ module derived_types
integer number_vis_eddies ! number of visible eddies
- integer, allocatable:: extr_map(:, :) ! (nlon, nlat) At a point
- ! of extremum SSH with out_cont found, of sufficient amplitude:
- ! identification number or this extremum. At a point of extremum
- ! SSH with no out_cont or out_cont of insufficient amplitude: the
- ! opposite of the identification number or this extremum. 0 at
- ! other points.
+ integer, allocatable:: extr_map(:, :) ! (nlon, nlat)
+ ! At a point of extremum SSH with out_cont found, of sufficient
+ ! amplitude: identification number or this extremum. At a point
+ ! of extremum SSH with no out_cont or out_cont of insufficient
+ ! amplitude: the opposite of the identification number or this
+ ! extremum. 0 at other points.
integer, allocatable:: ind_extr(:, :) ! (2, number_vis_eddies)
! indices in the global grid of ssh extrema
diff --git a/successive_overlap.f b/successive_overlap.f
new file mode 100644
index 00000000..78c1f44d
--- /dev/null
+++ b/successive_overlap.f
@@ -0,0 +1,81 @@
+module successive_overlap_m
+
+ implicit none
+
+contains
+
+ subroutine successive_overlap(unit_graph, nlon, nlat, flow, j, k)
+
+ ! Finds edges between flow(j - 1) and flow(j), corresponding to
+ ! dates k - 1 and k. Writes these edges to graph.txt. Updates
+ ! flow(j - 1)%list_vis%delta_out and flow(j)%list_vis%delta_in.
+
+ use contour_531, only: polyline
+ use derived_types, only: snapshot
+ use gpc_f, only: gpc_polygon_clip_f, GPC_INT, polygon
+ use spherical_polygon_area_m, only: spherical_polygon_area
+ use spherical_polyline_area_m, only: spherical_polyline_area
+ use weight_m, only: weight
+
+ integer, intent(in):: unit_graph ! logical unit for graph file
+ integer, intent(in):: nlon, nlat
+ type(snapshot), intent(inout):: flow(:) ! (max_delta + 1)
+ integer, intent(in):: j, k
+
+ ! Local:
+ integer i1, i1_lon, i1_lat, m, l, i2
+ type(polyline) polyline_1, polyline_2
+ type(polygon) res_pol
+
+ integer, parameter:: dist_lim = 12
+ ! We look for an overlapping eddy at dist_lim (in grid points) of
+ ! the first extremum.
+
+ !-----------------------------------------------------------------------
+
+ do i1 = 1, flow(j - 1)%number_vis_eddies
+ if (flow(j - 1)%list_vis(i1)%suff_amp) then
+ polyline_1 = merge(flow(j - 1)%list_vis(i1)%speed_cont%polyline, &
+ flow(j - 1)%list_vis(i1)%out_cont%polyline, &
+ flow(j - 1)%list_vis(i1)%speed_cont%n_points /= 0)
+ i1_lon = flow(j - 1)%ind_extr(1, i1)
+ i1_lat = flow(j - 1)%ind_extr(2, i1)
+
+ DO m = max(i1_lat - dist_lim, 1), min(i1_lat + dist_lim, nlat)
+ DO l = max(i1_lon - dist_lim, 1), min(i1_lon + dist_lim, nlon)
+ if (flow(j)%extr_map(l, m) > 0) then
+ i2 = flow(j)%extr_map(l, m)
+
+ if (flow(j)%list_vis(i2)%cyclone &
+ .eqv. flow(j - 1)%list_vis(i1)%cyclone) then
+ polyline_2 &
+ = merge(flow(j)%list_vis(i2)%speed_cont%polyline, &
+ flow(j)%list_vis(i2)%out_cont%polyline, &
+ flow(j)%list_vis(i2)%speed_cont%n_points /= 0)
+ call gpc_polygon_clip_f(GPC_INT, polygon(nparts = 1, &
+ part = [polyline_1], hole = [.false.]), &
+ polygon(nparts = 1, part = [polyline_2], &
+ hole = [.false.]), res_pol)
+
+ if (res_pol%nparts /= 0) then
+ ! polyline_1 and polyline_2 overlap
+ if (spherical_polygon_area(res_pol) >= 0.25 &
+ * min(spherical_polyline_area(polyline_1), &
+ spherical_polyline_area(polyline_2))) then
+ write(unit_graph, fmt = *) k - 1, i1, k, i2, &
+ weight(flow(j - 1)%list_vis(i1), &
+ flow(j)%list_vis(i2))
+ flow(j - 1)%list_vis(i1)%delta_out = 1
+ flow(j)%list_vis(i2)%delta_in = 1
+ end if
+ end if
+ end if
+ end if
+ end DO
+ end DO
+ end if
+ end do
+
+ end subroutine successive_overlap
+
+end module successive_overlap_m
--
GitLab