Commit bd370a6d authored by Robert Schweppe's avatar Robert Schweppe
Browse files

Merge branch 'add_functions_to_mo_poly' into 'main'

Add functions to mo_poly.f90 to modify any polygon to work with in_poly (from branch in chs/forces)

See merge request !55
parents c69786b7 91ba511b
Pipeline #99748 passed with stages
in 7 minutes and 11 seconds
......@@ -32,7 +32,7 @@ documentation:
stage: test
variables:
GIT_CLONE_PATH: $CI_BUILDS_DIR/$CI_RUNNER_SHORT_TOKEN/$CI_PROJECT_PATH/$CI_COMMIT_REF_NAME/$CI_JOB_NAME
GIT_FETCH_EXTRA_FLAGS: --all --no-tags --prune --quiet
GIT_FETCH_EXTRA_FLAGS: --no-tags --prune --quiet
script:
- source hpc-module-loads/eve.chs-conda01
- module load GCCcore/9.3.0 texlive/2020
......
......@@ -2,6 +2,15 @@
All notable changes to **FORCES** will be documented in this file.
## v0.4.0 - 2022-06
- See the git [diff](https://git.ufz.de/chs/forces/-/compare/v0.3.2...v0.4.0) for details.
### Enhancements
- mo_poly: added mo_poly.fypp replacement for mo_poly functions
- mo_poly: added new routines and tests
- `orientpoly` (calculate orientation of coords in polygon),
- `mod_pole` (modify coords of grid to include poles on Cartesian coord system) and
- `mod_shift` (shift longitude values by 180 degrees)
## v0.3.2 - 2022-06
See the git [diff](https://git.ufz.de/chs/forces/-/compare/v0.3.1...v0.3.2) for details.
......
......@@ -11,22 +11,23 @@ if(NOT FYPP)
message(WARNING "Preprocessor fypp not found!")
else()
message(STATUS "Preprocessor fypp found!")
# Create a list of the files to be preprocessed
set(fppFiles mo_netcdf.fypp)
# Create a list of the files to be preprocessed
set(fppFiles mo_netcdf.fypp mo_poly.fypp)
# Custom preprocessor flags
if(DEFINED CMAKE_MAXIMUM_RANK)
set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}")
elseif(f03rank)
set(fyppFlags)
# Fortran03 standard supports rank 15 arrays, activate here, if required
# elseif(f03rank)
# set(fyppFlags)
else()
set(fyppFlags "-DVERSION90")
endif()
# Pre-process: .fpp -> .f90 via Fypp
fypp_f90("${fyppFlags}" "${fppFiles}" fyppOutFiles)
list(FILTER sources EXCLUDE REGEX ".*mo_netcdf.f90$")
list(FILTER sources EXCLUDE REGEX ".*mo_poly.f90$")
endif()
# create library
add_library(${LIB_NAME} ${sources} ${fyppOutFiles})
target_include_directories(${LIB_NAME} PUBLIC ${CMAKE_CURRENT_BINARY_DIR})
......
This diff is collapsed.
!> \file mo_poly.f90
!> \copydoc mo_poly
#:include "common.fypp"
#:set KINDS_TYPES = REAL_KINDS_TYPES
#:set KINDS = REAL_KINDS
#:set PROCEDURES = ["areapoly", "center_of_mass", "inpoly", "orientpoly", "mod_pole", "mod_shift"]
!> \brief Polygon calculations.
!> \details
!! This module determines wether a 2D point lies inside, outside, or
!! on the vertice/edge of a 2D polygon
!! and is part of the UFZ CHS Fortran library.
!> \author Juliane Mai
!> \date Jul 2012
module mo_poly
! This module determines wether a 2D point lies inside, outside, or
! on the vertice/edge of a 2D polygon
! and is part of the UFZ CHS Fortran library.
!
! Written Juliane Mai, July 2012
! Modified Maren Goehler, July 2012 - area & center of mass
! License
! -------
! This file is part of the UFZ Fortran library.
! The UFZ Fortran library is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! The UFZ Fortran library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
! You should have received a copy of the GNU Lesser General Public License
! along with the UFZ Fortran library (LICENSE).
! If not, see <http://www.gnu.org/licenses/>.
! Copyright 2012 Juliane Mai
use mo_kind, only: i4, ${', '.join(KINDS)}$
use mo_utils, only: eq, ge, le, ne
implicit none
#:for proc in PROCEDURES
public :: ${proc}$
!> \copydoc mo_poly::${proc}$_${KINDS_TYPES[0][0]}$
interface ${proc}$
#:for kind, type in KINDS_TYPES
module procedure ${proc}$_${kind}$
#:endfor
end interface
#:endfor
! ------------------------------------------------------------------
private
! ------------------------------------------------------------------
contains
! ------------------------------------------------------------------
#:def areapoly_template(kind, type)
!> \brief Area of polygon
!> \details Function for computing the area of a polygon (2D, convex or not).
!! Function for computing the area of a polygon
!! The polygon can be convex or not.
!!
!! The method is only applicable for 2D polygons and points.
!!
!! \b Example
!!
!! \code{.f90}
!! polygon(:,1) = (/ 1.0_dp,2.0_dp,2.0_dp,1.0_dp /)
!! polygon(:,2) = (/ 1.0_dp,1.0_dp,2.0_dp,2.0_dp /)
!!
!! area = areapoly( polygon )
!!
!! ! --> area = 1
!! \endcode
!!
!! See also example in test directory.
!!
!! \b Literature
!!
!! 1. http://de.wikipedia.org/wiki/Geometrischer_Schwerpunkt
!!
!> \return Area of polygon
function areapoly_${kind}$(coord) result(areapoly)
!> coordinates of the polygon in question
${type}$, dimension(:,:), intent(in) :: coord
${type}$ :: areapoly
! local variables
integer(i4) :: i,k ! loop
integer(i4) :: nedges ! number of coordinates
${type}$ :: xsum ! for summing up
${type}$ :: ysum ! for summing up
xsum = 0.0_${kind}$
ysum = 0.0_${kind}$
nedges = size(coord,1)
do i = 1, nedges
if (i == nedges) then
k = 1_i4
else
k = i + 1_i4
end if
xsum = xsum + ( coord(i,1) * coord(k,2) )
ysum = ysum + ( coord(i,2) * coord(k,1) )
end do
areapoly = 0.5_${kind}$ * (xsum - ysum)
end function areapoly_${kind}$
#:enddef areapoly_template
#:def center_of_mass_template(kind, type)
!> \brief Center of mass of polygon.
!> \details Function for computing the center of mass of a polygon (2D, convex or not).
!!
!! \f[ A = \sum_{i}(x_{i}y_{i+1}-x_{i+1}y_{i}) \f]
!! \f[ x_s = \frac{1}{6A} \sum_i(x_i+x_{i+1})(x_iy_{i+1}-x_{i+1}y_i) \f]
!! \f[ y_s = \frac{1}{6A} \sum_i(y_i+y_{i+1})(x_iy_{i+1}-x_{i+1}y_i) \f]
!!
!! \b Example
!!
!! \code{.f90}
!! polygon(:,1) = (/ 1.0_dp,2.0_dp,2.0_dp,1.0_dp /)
!! polygon(:,2) = (/ 1.0_dp,1.0_dp,2.0_dp,2.0_dp /)
!!
!! com = center_of_mass( polygon )
!!
!! ! --> com = (/1.5_dp, 1.5_dp/)
!! \endcode
!!
!! See also example in test directory
!!
!! \b Literature
!!
!! 1. http://de.wikipedia.org/wiki/Geometrischer_Schwerpunkt
!!
!> \return Center of mass of polygon.
function center_of_mass_${kind}$(coord) result(center_of_mass)
!> coordinates of polygon in question
${type}$, dimension(:,:), intent(in) :: coord
${type}$, dimension(2) :: center_of_mass
! local variables
integer(i4) :: i,k ! loop
integer(i4) :: nedges ! number of coordinates
${type}$ :: area ! area of the polygon
${type}$ :: xsum ! for summing up
${type}$ :: ysum ! for summing up
xsum = 0.0_${kind}$
ysum = 0.0_${kind}$
nedges = size(coord,1)
area = areapoly_${kind}$(coord)
do i = 1, nedges
if (i == nedges ) then
k = 1_i4
else
k = i + 1_i4
end if
! multiply x coord by the y coord of next vertex
xsum = xsum + ((coord(i,1) + coord(k,1)) * &
((coord(i,1) * coord(k,2) - coord(k,1) * coord(i,2))))
ysum = ysum + ((coord(i,2) + coord(k,2)) * &
((coord(i,1) * coord(k,2) - coord(k,1) * coord(i,2))))
end do
center_of_mass(1) = 1.0_${kind}$ / (6.0_${kind}$ * area) * xsum
center_of_mass(2) = 1.0_${kind}$ / (6.0_${kind}$ * area) * ysum
end function center_of_mass_${kind}$
#:enddef center_of_mass_template
#:def inpoly_template(kind, type)
!> \brief Determination point of polygon.
!> \details Determines whether a 2D point is inside, outside or on vertex of a polygon (2D, convex or not).
!!
!! The original version of the source code (pnpoly) was implemented by
!! W. Randolph Franklin. It was insufficiently assigning vertex/edge points.
!!
!! \b Example
!!
!! \code{.f90}
!! polygon(:,1) = (/ 1.0_dp,2.0_dp,2.0_dp,1.0_dp /)
!! polygon(:,2) = (/ 1.0_dp,1.0_dp,2.0_dp,2.0_dp /)
!!
!! point = (/ 1.5, 1.5 /)
!!
!! call inpoly( point, polygon, inside )
!!
!! ! --> inside = 1 ... point is inside the polygon
!! \endcode
!!
!! See also example in test directory
!!
!! \b Literature
!!
!! 1. https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
!!
!> \return Whether point is inside (=1), outside (=-1) or on a vertex/edge of the polygon (=0)
subroutine inpoly_${kind}$(p,coord,erg)
!> point in question
${type}$, dimension(2), intent(in) :: p
!> coordinates of the polygon
${type}$, dimension(:, :), intent(in) :: coord
!> result:
!! inside: erg = 1
!! outside: erg = -1
!! on vertex/edge: erg = 0
integer(i4), intent(out) :: erg
! local variables
${type}$, dimension(size(coord,1)) :: x, y
${type}$ :: lx, ly
logical :: mx,my,nx,ny, test1, test2
integer(i4) :: n, i, j
n = size(coord,1)
do i=1,n
x(i)=coord(i,1)-p(1)
y(i)=coord(i,2)-p(2)
! check if point is equal to any coord
if ( eq(x(i),0.0_${kind}$) .and. eq(y(i),0.0_${kind}$) ) then
erg=0_i4
return
end if
end do
erg=-1_i4
do i=1,n
j=1+mod(i,n)
! vertical vertex
if ( eq(coord(i,1),coord(j,1)) .and. eq(coord(i,1),p(1)) ) then
ly = (p(2)-coord(j,2)) / (coord(i,2)-coord(j,2))
if ( ge(ly,0.0_${kind}$) .and. le(ly,1.0_${kind}$) ) then
erg=0_i4
return
end if
end if
! horizontal vertex
if ( eq(coord(i,2),coord(j,2)) .and. eq(coord(i,2),p(2)) ) then
lx = (p(1)-coord(j,1)) / (coord(i,1)-coord(j,1))
if ( ge(lx,0.0_${kind}$ ) .and. le(lx,1.0_${kind}$) ) then
erg=0_i4
return
end if
end if
!
mx = ge(x(i),0.0_${kind}$)
nx = ge(x(j),0.0_${kind}$)
my = ge(y(i),0.0_${kind}$)
ny = ge(y(j),0.0_${kind}$)
test1 = .not.((my.or.ny).and.(mx.or.nx)).or.(mx.and.nx)
test2 = .not.(my.and.ny.and.(mx.or.nx).and..not.(mx.and.nx))
if (.not. test1) then
if (test2) then
if ((y(i)*x(j)-x(i)*y(j))/(x(j)-x(i)) < 0.0_${kind}$) then
cycle
else
if ((y(i)*x(j)-x(i)*y(j))/(x(j)-x(i)) > 0.0_${kind}$) then
erg = -erg
cycle
else
erg = 0_i4
return
end if
end if
else
erg=-erg
end if
end if
end do
end subroutine inpoly_${kind}$
#:enddef inpoly_template
#:def orientpoly_template(kind, type)
!> \brief Check orientation of polygon
!> \details Function for checking the orientation of a polygon (2D, convex or not).
!> \return Integer indicating orientation (counter-clockwise: -1, clockwise: 1, none: 0)
function orientpoly_${kind}$(coord) result(orientpoly)
!> coordinates of the polygon in question
${type}$, dimension(:,:), intent(in) :: coord
integer(i4) :: orientpoly
!> result:
!! clockwise: orientpoly = 1
!! counter-clockwise: orientpoly = -1
!! flat line or similar irregular shape: orientpoly = 0
! local variables
integer(i4) :: n
${type}$ :: sum_edges
! calculate sum over the edges, (x2 - x1)(y2 + y1) as in
! https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order
n = size(coord, 1)
! use a vectorized version of sum over all (x2 -x1)*(y2-y1)
sum_edges = sum((coord(2:n, 1) - coord(1:n-1, 1)) * (coord(2:n, 2) + coord(1:n-1, 2)))
sum_edges = sum_edges + (coord(1, 1) - coord(n, 1)) * (coord(1, 2) + coord(n, 2))
if (eq(sum_edges, 0._${kind}$)) then
orientpoly = 0_i4
else if (sum_edges < 0._${kind}$) then
orientpoly = -1_i4
else
orientpoly = 1_i4
end if
end function orientpoly_${kind}$
#:enddef orientpoly_template
#:def mod_pole_template(kind, type)
!> \brief Modify polygon so it covers pole correctly
!> \details Modifies a polygon (2D, convex or not) to include pole when passed to inpoly
!! this function is intended to modify a given polygon, so it can be represented in a Cartesian grid
!! the use case is a polygon (e.g. 120,80, -120,80, 0,80 covering the north pole) that is not represented on
!! Cartesian lat-lon grid as polygon.
!! The script inserts additional coordinates, so the pole is covered (e.g. 180,80, 180,90, -180,90, -180,80)
!! See test cases for examples.
!> \return modified coordinates
function mod_pole_${kind}$(coord, meridian_arg) result(coord_mod)
!> coordinates of the polygon in question
${type}$, dimension(:,:), intent(in) :: coord
!> meridian that represents discontinuity, defaults to 180.0
${type}$, intent(in), optional :: meridian_arg
${type}$, dimension(:,:), allocatable :: coord_mod
! local variables
${type}$ :: meridian
${type}$ :: break, a
integer(i4) :: i, j, k, n
if (present(meridian_arg)) then
meridian = meridian_arg
else
meridian = 180._${kind}$
end if
n = size(coord, 1)
! determine location where meridian is crossed
! find the maximum and minimum longitudes
i = maxloc(coord(:, 1), 1)
j = minloc(coord(:, 1), 1)
! determine size of new coords array
k = n + 2
if (ne(coord(i, 1), meridian)) then
k = k + 1
end if
if (ne(coord(j, 1), meridian * (-1._${kind}$))) then
k = k + 1
end if
allocate(coord_mod(k, 2))
! polygon covers a pole
if (mod(i,n)+1 == j) then
! the coord pair after i is j, so longitudes are ascending, so north pole is contained
coord_mod(1:i, :) = coord(1:i, :)
k = i
! if the maxval is not meridian, we need to add point at intersection of meridian and points i and j
if (ne(coord(i, 1), meridian)) then
a = meridian - coord(i, 1)
break = coord(i, 2) + a / (a + abs(meridian + coord(j, 1))) * (coord(j, 2) - coord(i, 2))
coord_mod(k+1, :) = [meridian, break]
k = k + 1
end if
! add the points meridian,90 and meridian * -1, 90
coord_mod(k+1:k+2, 1) = [meridian, meridian * (-1._${kind}$)]
coord_mod(k+1:k+2, 2) = [90._${kind}$, 90._${kind}$]
k = k + 2
! if the minval is not meridian * -1, we need to add point at intersection of meridian and points i and j
if (ne(coord(j, 1), meridian * (-1._${kind}$))) then
a = meridian - coord(i, 1)
break = coord(i, 2) + a / (a + abs(meridian + coord(j, 1))) * (coord(j, 2) - coord(i, 2))
coord_mod(k+1, :) = [meridian * (-1._${kind}$), break]
k = k + 1
end if
! add the remaining coordinates
if (j > 1) coord_mod(k+1:k+1+n-j, :) = coord(j:n, :)
else if (mod(j,n)+1 == i) then
! the coord pair after j is i, so longitudes are descending, so south pole is contained
coord_mod(1:j, :) = coord(1:j, :)
k = j
! if the minval is not meridian * -1, we need to add point at intersection of meridian and points i and j
if (ne(coord(j, 1), meridian * (-1._${kind}$))) then
a = abs(meridian + coord(j, 1))
break = coord(j, 2) + a / (a + meridian - coord(i, 1)) * (coord(i, 2) - coord(j, 2))
coord_mod(k+1, :) = [meridian * (-1._${kind}$), break]
k = k + 1
end if
! add the points meridian * -1, -90 and meridian, -90
coord_mod(k+1:k+2, 1) = [meridian * (-1._${kind}$), meridian]
coord_mod(k+1:k+2, 2) = [-90._${kind}$, -90._${kind}$]
k = k + 2
! if the maxval is not meridian, we need to add point at intersection of meridian and points i and j
if (ne(coord(i, 1), meridian)) then
a = abs(meridian + coord(j, 1))
break = coord(j, 2) + a / (a + meridian - coord(i, 1)) * (coord(i, 2) - coord(j, 2))
coord_mod(k+1, :) = [meridian, break]
k = k + 1
end if
! add the remaining coordinates
if (i > 1) coord_mod(k+1:k+1+n-i, :) = coord(i:n, :)
! else: if there are multiple locations of minval or maxval, this edge case is not covered...
end if
end function mod_pole_${kind}$
#:enddef mod_pole_template
#:def mod_shift_template(kind, type)
!> \brief Shifts the (longitude) value 180 degrees
!> \details Modify a coordinate value
!> \return Shifted value
elemental function mod_shift_${kind}$(x_coord, meridian_arg) result(shifted)
!> coordinates of the polygon in question
${type}$, intent(in) :: x_coord
!> meridian that represents discontinuity, defaults to 180.0
${type}$, intent(in), optional :: meridian_arg
${type}$ :: shifted
${type}$ :: meridian
! shift values
if (present(meridian_arg)) then
meridian = meridian_arg
else
meridian = 180._${kind}$
end if
shifted = sign(abs(x_coord) - meridian, x_coord * (-1._${kind}$))
end function mod_shift_${kind}$
#:enddef mod_shift_template
#:for kind, type in KINDS_TYPES
$:areapoly_template(kind, type)
! ------------------------------------------------------------------
$:center_of_mass_template(kind, type)
! ------------------------------------------------------------------
$:inpoly_template(kind, type)
! ------------------------------------------------------------------
$:orientpoly_template(kind, type)
! ------------------------------------------------------------------
$:mod_pole_template(kind, type)
! ------------------------------------------------------------------
$:mod_shift_template(kind, type)
! ------------------------------------------------------------------
#:endfor
end module mo_poly
module test_mo_poly
use funit
use mo_kind, only: dp, sp, i4
use mo_poly, only: inpoly, areapoly, center_of_mass
use mo_message, only: error_message
use mo_kind, only: dp, i4
use mo_poly, only: inpoly, areapoly, center_of_mass, orientpoly, mod_pole, mod_shift
implicit none
integer(i4) :: N=4
real(dp), dimension(:,:), allocatable :: coord_dp
real(sp), dimension(:,:), allocatable :: coord_sp
real(dp), dimension(4, 2) :: coord_dp = reshape([&
1.0_dp,2.0_dp,2.0_dp,1.0_dp, &
1.0_dp,1.0_dp,2.0_dp,2.0_dp], &
shape(coord_dp))
integer(i4) :: inside
real(dp) :: area_dp
real(sp) :: area_sp
real(dp),dimension(2) :: com_dp
real(sp),dimension(2) :: com_sp
contains
! Test coordinate of a point compared to polygon
@test
subroutine test_poly_inpoly_dp()
allocate (coord_dp(N,2))
coord_dp(:,1) = (/ 1.0_dp,2.0_dp,2.0_dp,1.0_dp /)
coord_dp(:,2) = (/ 1.0_dp,1.0_dp,2.0_dp,2.0_dp /)
call inpoly( (/1.5_dp,1.5_dp/) , coord_dp, inside)
@assertEqual(inside, 1_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)
print*, 'THE POINT IS INSIDE THE POLYGON'
case(0)
print*, 'THE POINT IS ON AN EDGE OR AT A VERTEX'
end select
call inpoly( (/0.5_dp,1.5_dp/) , coord_dp, inside)
@assertEqual(inside, -1_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)
print*, 'THE POINT IS INSIDE THE POLYGON'
case(0)
print*, 'THE POINT IS ON AN EDGE OR AT A VERTEX'
end select
call inpoly( (/1.5_dp,1.0_dp/) , coord_dp, inside)
@assertEqual(inside, 0_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)
print*, 'THE POINT IS INSIDE THE POLYGON'
case(0)
print*, 'THE POINT IS ON AN EDGE OR AT A VERTEX'
end select
call inpoly( (/1.0_dp,1.5_dp/) , coord_dp, inside)
@assertEqual(inside, 0_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)
print*, 'THE POINT IS INSIDE THE POLYGON'
case(0)
print*, 'THE POINT IS ON AN EDGE OR AT A VERTEX'
end select
call inpoly( (/1.0_dp,1.0_dp/) , coord_dp, inside)
@assertEqual(inside, 0_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)
print*, 'THE POINT IS INSIDE THE POLYGON'
case(0)
print*, 'THE POINT IS ON AN EDGE OR AT A VERTEX'
end select
end subroutine test_poly_inpoly_dp
@test
subroutine test_poly_inpoly_sp()
allocate (coord_sp(N,2))
coord_sp(:,1) = (/ 1.0_sp,2.0_sp,2.0_sp,1.0_sp /)
coord_sp(:,2) = (/ 1.0_sp,1.0_sp,2.0_sp,2.0_sp /)
call inpoly( (/1.5_sp,1.5_sp/) , coord_sp, inside)
@assertEqual(inside, 1_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)
print*, 'THE POINT IS INSIDE THE POLYGON'
case(0)
print*, 'THE POINT IS ON AN EDGE OR AT A VERTEX'
end select
call inpoly( (/0.5_sp,1.5_sp/) , coord_sp, inside)
@assertEqual(inside, -1_i4)
select case (inside)
case(-1)
print*, 'THE POINT IS OUTSIDE THE POLYGON'
case(1)