https://github.com/cran/fields
Tip revision: 328f151e07864e976578f9134025cbfdb72eb2f3 authored by Doug Nychka on 19 July 2013, 23:45:29 UTC
version 6.8
version 6.8
Tip revision: 328f151
igpoly.f
subroutine igpoly(nx, xg,ny,yg,np,xp,yp,ind)
!----------------------------------------------------------------------
! This subroutine determines whether or not an 2-d point (xg(i),yg(j))
! element is inside a (closed) set of points xp,yp that
! are assumed to form polygon.
! result is an matrix indicating comparision for all
! combinations of xg and yg
!----------------------------------------------------------------------
integer np ! # of points in polygon
integer nx,ny ! # points to check
real xg(nx) ! x grid values to check
real yg(ny) ! y grid values to check
real xp(np) ! 2d-locations of polygon
real yp(np)
real x1, x2, y1,y2 ! min and max of x and y
real temp, xt, yt
integer ind(nx,ny) ! THE ANSWER : ind(i)=1 if point xd(i),yd(i) is
! in polygon 0 otherwise
integer in
x1= xp(1)
x2= xp(2)
y1= yp(1)
y2= yp(1)
!
! find the minima and maxima of the polygon coordinates
! i.e. the smallest rectangle containing the polygon.
do j = 1,np
temp= xp(j)
if( temp.lt.x1) then
x1 = temp
endif
if( temp.gt.x2) then
x2 = temp
endif
temp= yp(j)
if( temp.lt.y1) then
y1 = temp
endif
if( temp.gt.y2) then
y2 = temp
endif
enddo
! loop over all
! grid points
do i = 1, nx
do j = 1,ny
xt= xg(i)
yt= yg(j)
! quick test that point is inside the bounding rectangle
! of the polygon
if( (xt.le. x2).and. (xt.ge.x1).and.
* (yt.le.y2).and.(yt.ge.y1) ) then
call inpoly2(xt,yt,np,xp,yp, in)
ind(i,j)=in
else
ind(i,j)= 0
endif
enddo
enddo
return
end