points-in-polygon routine (splus)

1. (get routines from Mark Durst's directory)

   rlogin cedr -l merrill
   cd $MDOCS/splus/ptsinpoly
   cp -p /h/durst/progs/sprogs/dwm.polyg/ptsinpoly.s .
   more ptsinpoly.s
   cp -p /h/durst/progs/sprogs/dwm.polyg/ptsinpoly.f .
   more ptsinpoly.f
   ls -alr /h/durst/progs/sprogs/dwm.polyg/
   

   Mark Durst's directory contains:

   /h/durst/progs/sprogs/dwm.polyg/:
   total 262
   drwxr-sr-x  3 durst         512 Apr 12 08:46 .
   drwxr-sr-x  5 durst         512 Apr 11 14:39 ..
   drwxr-sr-x  2 durst         512 Apr 12 11:03 .Data
   -rwxr-xr-x  1 durst      163840 Apr 12 08:46 a.out
   -rwxr-xr-x  1 durst        4527 Apr 12 08:46 mytest.out
   -rwxr-xr-x  1 durst        6315 Apr 12 08:46 pt_test.f
   -rwxr-xr-x  1 durst        6320 Apr 11 16:38 pt_test.f.0
   -rwxr-xr-x  1 durst        6321 Apr 12 08:31 pt_test.f.1
   -rwxr-xr-x  1 durst        6321 Apr 11 16:40 pt_test.f.~1~
   -rwxr-xr-x  1 durst         278 Apr 12 08:20 ptinpoly.s
   -rwxr-xr-x  1 durst        4310 Apr 12 07:58 ptsinpoly.f
   -rwxr-xr-x  1 durst        7284 Apr 11 15:49 ptsinpoly.f.0
   -rwxr-xr-x  1 durst        4075 Apr 12 07:42 ptsinpoly.f.bak
   -rwxr-xr-x  1 durst        2100 Apr 12 07:58 ptsinpoly.o
   -rwxr-xr-x  1 durst         332 Apr 12 07:37 ptsinpoly.s
   -rwxr-xr-x  1 durst       13500 Apr 11 16:42 test.out
   -rwxr-xr-x  1 durst       22527 Apr 11 17:00 test2.out
   
   /h/durst/progs/sprogs/dwm.polyg/.Data:
   total 76
   drwxr-sr-x  2 durst         512 Apr 12 11:03 .
   drwxr-sr-x  3 durst         512 Apr 12 08:46 ..
   -rwxr-xr-x  1 durst       37437 Apr 12 11:03 .Audit
   -rwxr-xr-x  1 durst          20 Apr 12 11:03 .Last.value
   -rwxr-xr-x  1 durst          64 Apr 12 08:03 .Random.seed
   -rwxr-xr-x  1 durst         285 Apr 12 08:58 last.dump
   -rwxr-xr-x  1 durst         933 Apr 12 08:20 ptinpoly
   -rwxr-xr-x  1 durst        1148 Apr 12 07:31 ptsinpoly
   -rwxr-xr-x  1 durst        1616 Apr 12 08:00 tmp
   -rwxr-xr-x  1 durst       20063 Apr 12 08:08 tmp.mat
   -rwxr-xr-x  1 durst         592 Apr 12 09:53 tmp.out
   -rwxr-xr-x  1 durst        1168 Apr 12 09:49 tmp.x
   -rwxr-xr-x  1 durst          80 Apr 12 08:20 tmp.xy
   -rwxr-xr-x  1 durst        1168 Apr 12 09:50 tmp.y
   -rwxr-xr-x  1 durst         416 Apr 12 10:01 tmp2.out
   -rwxr-xr-x  1 durst         816 Apr 12 09:59 tmp2.x
   -rwxr-xr-x  1 durst         816 Apr 12 10:00 tmp2.y
   

2. compile and load source:

   rlogin parep2 -l merrill
   cd $MDOCS/splus/ptsinpoly
   module load lang
   f77 -c ptsinpoly.f
   splus
   source ("ptsinpoly.s")
   ptsinpoly
   

   Result is:

   ptsinpoly
   function(x, y, xinf, yinf, xy)
   {
           z <- .Fortran("ptsinpoly",
                   n.check = as.integer(length(x)),
                   as.single(x),
                   as.single(y),
                   as.single(xinf),
                   as.single(yinf),
                   n.vert = as.integer(length(xy)/2),
                   as.single(xy),
                   iotmsg = as.integer(0),
                   inpoly = as.integer(rep(0, length(x))))
           z$inpoly
   }
   

3. execute:

   rlogin parep2 -l merrill
   cd $MDOCS/splus/ptsinpoly
   splus
   dyn.load("ptsinpoly.o")
   

4. square with corners at +/- 1; (xinf,yinf) lies outside:

   xinf_100
   yinf_100
   xy_c(-1,-1, -1,1, 1,1, 1,-1)
   xinf
   yinf
   xy
   

   result is:

   >xinf
   [1] 100
   >yinf
   [1] 100
   >xy
   [1] -1 -1 -1  1  1  1  1 -1
   

5. three points, two in and one out:

   x_c(.5,.8,1.2)
   y_c(.4,.7,1.3)
   inpoly_ptsinpoly(x,y,xinf,yinf,xy)
   x
   y
   inpoly
   

   result is:

   >x
   [1] 0.5 0.8 1.2
   >y
   [1] 0.4 0.7 1.3
   >inpoly
   [1] 1 1 0
   

6. 10 random points in a square with corners at +/- 2:

   npts_10
   a_rep(1,npts)
   x_-2+2*runif(a)
   y_-2+2*runif(a)
   inpoly_ptsinpoly(x,y,xinf,yinf,xy)
   x
   y
   inpoly
   npts
   sum(inpoly)
   

   result is:

   > x
    [1] -0.6871102 -1.9059346 -0.4147711 -1.4099114 -0.2549921 -1.1717878
    [7] -1.2778158 -0.8242569 -1.6333318 -0.3460031
   > y
    [1] -0.2250555 -1.1036882 -0.9963645 -0.6585269 -1.3451567 -1.1078135
    [7] -0.9130003 -0.8798343 -1.8819547 -1.4788028
   > inpoly
    [1] 1 0 1 0 0 0 0 1 0 0
   > npts
   [1] 10
   > sum(inpoly)
   [1] 3
   

7. 200000 random points in a square with corners at +/- 2:

   npts_200000
   a_rep(1,npts)
   x_-2+2*runif(a)
   y_-2+2*runif(a)
   inpoly_ptsinpoly(x,y,xinf,yinf,xy)
   npts
   sum(inpoly)
   

   result is:

   					time on cedr	time on parep2
   > a_rep(1,npts) 			8 secs		2 secs
   > x_-2+2*runif(a)			20 secs		4 secs
   > y_-2+2*runif(a)			20 secs		4 secs
   > inpoly_ptsinpoly(x,y,xinf,yinf,xy)	40 secs		15 secs
   > npts
   [1] 2e+05
   > sum(inpoly)
   [1] 50094