Orient the polygon rings

Common/CMakeLists.txt

 add_subdirectory(Tests)
 target_sources(test_get_1_outerm PRIVATE derived_types.f90
-  spher_polyline_area.f90 shpc_create.f90 write_eddy.f90 shpc_close.f90)
+  spher_polyline_area.f90 shpc_create.f90 write_eddy.f90 shpc_close.f90
+  ccw_orient.f90)
 target_sources(test_set_all_outerm PRIVATE derived_types.f90
-  spher_polyline_area.f90 shpc_create.f90 write_eddy.f90 shpc_close.f90)
+  spher_polyline_area.f90 shpc_create.f90 write_eddy.f90 shpc_close.f90
+  ccw_orient.f90)
 target_sources(test_max_speed_contour_ssh PRIVATE derived_types.f90)
 target_sources(test_nearby_extr PRIVATE derived_types.f90 read_snapshot.f90
   read_eddy.f90 read_field_indices.f90 shpc_open.f90 shpc_close.f90)
 target_sources(test_set_max_speed PRIVATE derived_types.f90
   spher_polyline_area.f90 shpc_open.f90 shpc_close.f90 read_field_indices.f90
-  read_eddy.f90 shpc_create.f90 write_eddy.f90)
+  read_eddy.f90 shpc_create.f90 write_eddy.f90 ccw_orient.f90)
 target_sources(test_spher_polyline_area PRIVATE spher_polyline_area.f90)
 target_sources(inst_eddies PRIVATE write_eddy.f90 spher_polyline_area.f90
   derived_types.f90 shpc_create.f90 shpc_close.f90 shpc_open.f90
-  read_field_indices.f90)
+  read_field_indices.f90 ccw_orient.f90)
 target_sources(test_write_eddy PRIVATE derived_types.f90 shpc_open.f90
   shpc_close.f90 read_field_indices.f90 shpc_create.f90)
 

Common/ccw_orient.f90

+module ccw_orient_m
+
+  implicit none
+
+contains
+
+  subroutine ccw_orient(c)
+
+    ! Revert orientation of contour if area is negative. We are
+    ! choosing the GeoJSon convention: counterclockwise external ring,
+    ! over the ESRI shapefile convention, which is opposite. The main
+    ! point is probably to follow a convention rather than let the
+    ! polygon have random order.
+
+    use derived_types, only: ssh_contour
+
+    type(ssh_contour), intent(inout):: c
+
+    !---------------------------------------------------------------------
+
+    if (c%area < 0) then
+       ! Reverse orientation:
+       c%polyline%points = c%polyline%points(:, c%polyline%n_points:1:- 1)
+       c%area = - c%area
+    end if
+
+  end subroutine ccw_orient
+
+end module ccw_orient_m

Common/spher_polyline_area.f90

@@ -7,8 +7,9 @@ contains
   pure real function spher_polyline_area(p)
 
     ! 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.
+    ! area (positive if vertices are in counter-clockwise order) of a
+    ! polyline in longitude, sin(latitude) with the same
+    ! vertices. Result in m2.
 
     ! Libraries:
     use contour_531, only: polyline
@@ -27,7 +28,7 @@ contains
 
     v(1, :) = p%points(1, :)
     v(2, :) = sin(p%points(2, :))
-    spher_polyline_area = r_Earth**2 * abs(polygon_area_2d(v))
+    spher_polyline_area = r_Earth**2 * polygon_area_2d(v)
 
   end function spher_polyline_area
 

Inst_eddies/get_1_outerm.f90

@@ -26,6 +26,7 @@ contains
     use contour_531, only: polyline
     use jumble, only: assert
 
+    use ccw_orient_m, only: ccw_orient
     use derived_types, only: ssh_contour, missing_ssh
     use good_contour_m, only: good_contour
     use spher_polyline_area_m, only: spher_polyline_area
@@ -105,6 +106,7 @@ contains
        end if
 
        get_1_outerm%area = spher_polyline_area(get_1_outerm%polyline)
+       call ccw_orient(get_1_outerm)
     end if
 
   end function get_1_outerm

Inst_eddies/set_max_speed.f90

@@ -28,6 +28,7 @@ contains
     use contour_531, only: convert_to_ind
     use geometry, only: polygon_point_dist_2d
 
+    use ccw_orient_m, only: ccw_orient
     use derived_types, only: eddy, null_ssh_contour, missing_ssh
     use good_contour_m, only: good_contour
     use inside_4_m, only: inside_4
@@ -107,6 +108,7 @@ contains
                 ! This may happen when the outermost contour is near land.
                 ! Stick to the contour coming from max_speed_contour_ssh.
                 e%speed_cont%area = spher_polyline_area(e%speed_cont%polyline)
+                call ccw_orient(e%speed_cont)
              else
                 ! Note the following test should raise an invalid exception
                 ! if e%max_speed is NaN.
@@ -118,6 +120,7 @@ contains
                    ! Stick to the contour coming from max_speed_contour_ssh.
                    e%speed_cont%area &
                         = spher_polyline_area(e%speed_cont%polyline)
+                   call ccw_orient(e%speed_cont)
                 end if
              end if
           end if

Overlap/overlap.f90

@@ -112,8 +112,8 @@ contains
              if (res_pol%nparts /= 0) then
                 ! polyline_1 and polyline_2 overlap
                 if (spher_polygon_area(res_pol) >= 0.25 &
-                     * min(spher_polyline_area(polyline_1), &
-                     spher_polyline_area(polyline_2))) then
+                     * min(abs(spher_polyline_area(polyline_1)), &
+                     abs(spher_polyline_area(polyline_2)))) then
                    write(unit_edge, fmt = *) (k - delta) * e_overestim + i1, &
                         k * e_overestim + i2
                    flow(j - delta)%list(i1)%delta_out &

Overlap/spher_polygon_area.f90

@@ -24,7 +24,7 @@ contains
 
     !------------------------------------------------------
 
-    forall(i = 1:p%nparts) parts_area(i) = spher_polyline_area(p%part(i))
+    forall(i = 1:p%nparts) parts_area(i) = abs(spher_polyline_area(p%part(i)))
     spher_polygon_area = sum(merge(- parts_area, parts_area, p%hole))
 
   end function spher_polygon_area