test_spher_polygon_area.f90

program test_spher_polygon_area

  ! Libraries:
  use gpc_f, only: polygon, shpobj2pol
  use jumble, only: get_command_arg_dyn
  use nr_util, only: deg_to_rad
  use shapelib, only: shpfileobject, shpobject, shpclose, shpdestroyobject
  use shapelib_03, only: shp_open_03, shp_read_object_03, shp_get_info_03
  
  use spher_polygon_area_m, only: spher_polygon_area

  implicit none

  character(len = :), allocatable:: filename
  TYPE(shpfileobject) hshp
  TYPE(shpobject) psobject
  type(POLYGON), allocatable:: p(:)
  type(POLYGON) merged_p
  integer ishape, j, nentities

  !---------------------------------------------------------------------

  call get_command_arg_dyn(1, filename, "Required argument: shapefile")
  call shp_open_03(hshp, filename, "rb")
  call shp_get_info_03(hshp, nentities)
  allocate(p(0:nentities - 1))
  
  do ishape = 0, nentities - 1
     call shp_read_object_03(hshp, ishape, psobject)
     p(ishape) = shpobj2pol(psobject)
     call shpdestroyobject(psobject)
  end do

  call shpclose(hshp)
  merged_p%nparts = sum(p%nparts)
  allocate(merged_p%part(merged_p%nparts), merged_p%hole(merged_p%nparts))

  j = 1
  do ishape = 0, nentities - 1
     merged_p%part(j:j + p(ishape)%nparts - 1) = p(ishape)%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(ishape)%nparts - 1) = .true.
     
     j = j + p(ishape)%nparts
  end do

  forall (j = 1:merged_p%nparts) &
       merged_p%part(j)%points = merged_p%part(j)%points * deg_to_rad
  print *, "Area = ", spher_polygon_area(merged_p) / 1e6, "km2"

end program test_spher_polygon_area