program findlonelypts
c - Cannibalized from "findoutliers", does the following:
c - Given:
c - A "vector file" (preferably, and this program is built specifically for this, of "RZ" variety)
c - A radius
c - A threshold, in meters
c - Search the entire vector file, and spit out all points with NO neighbors inside the radius that exceed the threshhold
c - VERTCON 3:
c LON LAT Meters PID
c (rvalm)
c 278.2166665250 24.5500001417 -0.415 AA0001
c 278.1999999639 24.5500001111 -0.408 AA0002
c 278.1946555556 24.5606000000 -0.408 AA0003
c 278.1938055556 24.5600000000 -0.409 AA0004
implicit real*8(a-h,o-z)
parameter(max = 500000)
character*80 card0,card(max)
character*200 fname,c
real*8 glo(max),gla(max)
real*8 qz(max),qz2(max)
real*8 sval(max)
real*8 rvals(max),rvals2(max)
real*8 rvalm(max),rvalm2(max)
character*6 xpid,pid(max),pidskip(max),pidhold(max)
logical*1 ldrop(max)
character*10 region,olddtm,newdtm
character*3 typ
logical allnoaz,noaz,badmag,badaz,swe
character*1 yesno,outform
integer*4 ilo(10),ihi(10)
logical*1 front
write(6,1)
1 format('Program findlonelypts')
write(6,2)
2 format('Enter file name to be searched: ',$)
read(5,'(a)')fname
c----------------
c - From the file name, get the region, and the old and new datums, as well as type
c - of data (lat, lon, hor or eht)
c----------------
idot = 0
ilo(1) = 1
do 3 i=1,len(trim(fname))
if(fname(i:i).eq.'.')then
idot = idot + 1
ihi(idot) = i-1
ilo(idot+1) = i+1
endif
3 continue
ihi(idot+1) = len(trim(fname))
typ = 'oht'
olddtm = fname(ilo(2) :ihi(2))
newdtm = fname(ilo(3) :ihi(3))
region = fname(ilo(4) :ihi(4))
write(6,4)
4 format('Enter search distance in arcMINUTES : ',$)
read(5,*)dxmin
c - convert to arcdegrees:
dx = dxmin / 60.d0
allnoaz = .true.
c - New for VERTCON3, dealing with (sometimes) residual stuff:
write(6,107)
107 format('Print output only if the actual ',
*'absolute value exceeds some threshhold? : ',$)
read(5,'(a)')yesno
if(yesno.eq.'y' .or. yesno.eq.'Y')then
write(6,108)
108 format(' Enter minimum absolute value to exceed: ',$)
read(5,*)xminabs
else
xminabs = 0
endif
98 write(6,97)
97 format('Search workedits file for newly ',
*'identified outliers to ignore? : ',$)
read(5,'(a)')yesno
if(yesno.ne.'y' .and. yesno.ne.'Y' .and.
*yesno.ne.'n' .and. yesno.ne.'N')goto 98
swe = .false.
if(yesno.eq.'y' .or. yesno.eq.'Y')swe = .true.
c - 2016/05/13
140 write(6,139)
139 format('Output format: standard(s) or workedits(w)? : ',$)
read(5,'(a)')yesno
if(yesno.eq.'s' .or. yesno.eq.'S')then
outform = 's'
elseif(yesno.eq.'w' .or. yesno.eq.'W')then
outform = 'w'
else
goto 140
endif
c----------------
c - Open the vector file, and report
c----------------
write(6,7)trim(fname)
7 format('Opening file: ',a)
open(1,file=fname,status='old',form='formatted')
write(6,6)trim(typ),trim(olddtm),trim(newdtm),trim(region)
6 format(
*'Type : ',a,/,
*'OldDtm : ',a,/,
*'NewDtm : ',a,/,
*'Region : ',a)
c----------------
c - Using the olddtm,newdtm and region, read in
c - any "workedits" points that've been changed since
c - we last created the vector file, to make sure
c - to ignore them.
c----------------
if(swe)then
write(6,5)
5 format('Searching workedits for additional PIDs to skip...')
open(90,file='Work/workedits',status='old',form='formatted')
c - Read each edit card into RAM, without regard for its content
701 read(90,702,end=703)c
if(len(c).eq.0)goto 701
if(c(1:1).eq.'#' .or. c(1:1).eq.' ')goto 701
c - Match on olddtm or get out
if(trim(c( 1: 10)).ne.trim(olddtm))goto 701
c - Match on newdtm or get out
if(trim(c( 12: 21)).ne.trim(newdtm))goto 701
c - Match on region or get out
if(trim(c( 23: 32)).ne.trim(region))goto 701
c - At this point, I have a match
neditsTotal = neditsTotal + 1
pidskip(neditsTotal) = c( 34: 39)
goto 701
702 format(a)
703 continue
c write(6,704)neditsTotal
704 format('Total Manual Edits Found: ',i6)
endif
101 format(f16.10,1x,f15.10,1x,6x ,1x,
*12x ,1x,9x ,1x,f9.3,1x,a6)
103 format('Bad Vector Length ',
*f16.10,1x,f15.10,1x,f6.2,1x,
*f12.2,1x,f9.5,1x,f9.3,1x,a6)
c----------------------------------
c Read entire vector file into RAM
c----------------------------------
i = 0
100 read(1,'(a)',end=102)card0
c write(6,'(a)')card0
c - Before counting this record as good, check against the "workedits" file, if we were
c - told to do so, and make sure there isn't a NEW outlier in that file that we haven't
c - yet swept up into the vector file through the "doit.bat" and "doit2.bat" procedures.
if(swe)then
do 231 iedits = 1,neditsTotal
if(card0(74:79).eq.pidskip(iedits))then
write(6,232)pidskip(iedits)
goto 100
endif
231 continue
endif
232 format('New workedits entry found, skipped: ',a6)
i = i + 1
card(i) = card0
read(card0,101)glo(i),gla(i),rvalm(i),pid(i)
rvalm2(i)= rvalm(i)**2
ldrop(i) = .false.
goto 100
102 npts = i
write(6,8)npts
8 format('Done reading file. ',/,
*6x,'Total pts found : ',i10)
c - Header of standard output
if(outform.eq.'s')write(6,9)
c----------------------------------
c Brute force: Each point against
c its neighbors. When searching
c for outliers, compare in meters.
c----------------------------------
k1 = 0
k2 = 0
do 300 i=1,npts
c write(6,*) ' Analyzing point # ',i
aveval = 0.d0
rmsval = 0.d0
noaz = .false.
if(allnoaz)noaz = .true.
badmag = .false.
badaz = .false.
c write(6,*) ' i = ',i
do 301 j=1,npts
c write(6,*) ' j = ',j
c - If this is the same point, skip
if(i.eq.j)goto 301
c - If I find a SINGLE point INSIDE the search radius, skip out of here and
c - begin searching the next "i-th" point
dlo = dabs(glo(j) - glo(i))
dla = dabs(gla(j) - gla(i))
if(dlo.lt.dx .and. dla.lt.dx)goto 300
301 continue
c - If I get here, this point has NO neighbors inside the search radius. Now
c - let's see if it exceeds our threshold.
k1 = k1 + 1
c write(6,401)pid(i)
c 401 format(' No neighbors! : ',a)
if(dabs(rvalm(i)) .lt. xminabs)goto 300
c - If I get here, report this as a lonely point that exceeded our threshold:
k2 = k2 + 1
if(outform.eq.'s')then
write(6,918)dxmin,xminabs,k2,i,gla(i),glo(i),
* rvalm(i),pid(i)
else
write(6,926)
* adjustl(olddtm),adjustl(newdtm),
* adjustl(region),pid(i),
* dxmin,xminabs,rvalm(i)
endif
read(card0,101)glo(i),gla(i),rvalm(i),pid(i)
300 continue
9 format(63x,
*'Count',12x,'Pt #',7x,'Lat',11x,'Lon',14x,'Val',4x,'PID')
918 format(' Lonely Point (RadArcMin: ',f5.1,' , ',
*' MinVal: ',f6.2,') -- ',i15,1x,i15,1x,f14.10,1x,f14.10,
*1x,f10.3,1x,a6)
926 format(3(a10,'|'),a6,'|',
* 'findlonelypts: RadMin=',i5.1,' and MinVal = ',f5.2,
* ' Val = ',f8.2)
write(6,900)trim(fname),dxmin,xminabs,npts,k1,k2
900 format('File: ',a,1x,'Radius,Arcmin: ',f5.1,1x,
*'MinAbs Value: ',f5.2,/,
*'Pts Found : ',i10,/,
*' Of these, # without neighbors : ',i10,/,
*' Of these, # over threshhold : ',i10)
end