In program test_spherical_polygon_area, allow reading a

multi-polygon. This is a useful test because, in successive_overlap,
the intersection of two polygons may give a multi-polygon.

Tests/test_spherical_polygon_area.f

 program test_spherical_polygon_area
 
+  use, intrinsic:: ISO_C_BINDING
+
   use gpc_f, only: polygon, shpobj2pol
   use jumble, only: get_command_arg_dyn
-  use nr_util, only: pi, assert
-  use shapelib, only: shpfileobject, shpobject, shpclose, shpdestroyobject
+  use nr_util, only: pi
+  use shapelib, only: shpfileobject, shpobject, shpclose, shpdestroyobject, &
+       shpgetinfo
   use shapelib_03, only: shp_open_03, shp_read_object_03
   use spherical_polygon_area_m, only: spherical_polygon_area
 
@@ -12,22 +15,44 @@ program test_spherical_polygon_area
   character(len = :), allocatable:: filename
   TYPE(shpfileobject) hshp
   TYPE(shpobject) psobject
-  type(POLYGON) p
+  type(POLYGON), allocatable:: p(:)
+  type(POLYGON) merged_p
   real, parameter:: deg2rad = pi / 180.
-  integer i
+  integer i, j, nentities, shapetype, dbffieldcount, dbfrecordcount
+  real(c_double) minbound(4), maxbound(4)
 
   !---------------------------------------------------------------------
 
   call get_command_arg_dyn(1, filename, "Required argument: shapefile")
   call shp_open_03(filename, "rb", hshp)
-  call shp_read_object_03(hshp, 0, psobject)
+  call shpgetinfo(hshp, nentities, shapetype, minbound, maxbound, &
+       dbffieldcount, dbfrecordcount)
+  allocate(p(nentities))
+  
+  do i = 1, nentities
+     call shp_read_object_03(hshp, i - 1, psobject)
+     p(i) = shpobj2pol(psobject)
+     call shpdestroyobject(psobject)
+  end do
+
   call shpclose(hshp)
-  call assert(psobject%nparts >= 1, "psobject%nparts >= 1")
-  call assert(psobject%panpartstart(1) == 0, "psobject%panpartstart(1)")
-  p = shpobj2pol(psobject)
-  call shpdestroyobject(psobject)
+  merged_p%nparts = sum(p%nparts)
+  allocate(merged_p%part(merged_p%nparts), merged_p%hole(merged_p%nparts))
+
+  j = 1
+  do i = 1, nentities
+     merged_p%part(j:j + p(i)%nparts - 1) = p(i)%part
+
+     ! Assume that the first part of each shape is the exterior ring,
+     ! and the following parts are holes (GeoJSON convention):
+     merged_p%hole(j) = .false.
+     merged_p%hole(j + 1:j + p(i)%nparts - 1) = .true.
+     
+     j = j + p(i)%nparts
+  end do
 
-  forall (i = 1:p%nparts) p%part(i)%points = p%part(i)%points * deg2rad
-  print *, "Area = ", spherical_polygon_area(p) / 1e6, "km2"
+  forall (i = 1:merged_p%nparts) &
+       merged_p%part(i)%points = merged_p%part(i)%points * deg2rad
+  print *, "Area = ", spherical_polygon_area(merged_p) / 1e6, "km2"
 
 end program test_spherical_polygon_area

Tests/tests.json

@@ -198,5 +198,11 @@
 	"args": ["$compil_prod_dir/test_spherical_polygon_area",
 		 "$input_dir/triangle"],
 	"title": "Spherical_polygon_area"
+    },
+    {
+	"args": ["$compil_prod_dir/test_spherical_polygon_area",
+		 "$input_dir/tri_square_hole"],
+	"title": "Area_multi_polygon",
+	"description": "Area of a multipolygon, with a hole in one polygon."
     }
 ]