program sift
c
c accesses a data bank of boundary lines, faults, etc. to get a 'lines file' for a particular area of interest.
c also winnows the number of line segments so that the detail preserved is appropriate to the scale.
c written by David Oppenheimer 2/86
c revised 5/12/93 to include nevada faults, david castillo line files, and reorganize code
c revised 7/15/94 to include 92 Jennings data, variable output format, ability to open new output file
c revised 3/08/95 to output MIF/MID format for Map Info
c
c Although this program has been used by the USGS, no warranty, expressed or implied, is made by the
c USGS or the United States Government as to the accuracy and functioning of the program and related
c program material not shall the fact of distribution constitute any such warranty, and no responsibility
c is assumed by the USGS in connection therewith.
c
integer ncode
c ! dimension of filname, coddat
parameter (ncode = 35)
c
integer ounit
c ! logical unit number of terminal output file
c ------------------------------------------------------------------------------------------
c THE FOLLOWING LINES NEED TO BE CHANGED TO REFLECT YOUR PARTICULAR SITE CHARACTERISTICS
c ------------------------------------------------------------------------------------------
character*22 prefix1
c ! pathname prefix
character*22 prefix2
c ! special path to handle copyrighted CDMG files
parameter (ounit = 0)
c ! UNIX version
c parameter (ounit = 6)
c ! VMS version
parameter (prefix1 = '/ftp/pub/outgoing/map/')
parameter (prefix2 = '/we/geysers/oppen/map/')
c ! UNIX version
c parameter (prefix = 'gsvax0::arc:[oppenheimer.map]')
c ! VMS version
c
character*1 ans
c ! answer
real azm
c ! azimuth
real bazm
c ! back azimuth
integer bunit
c ! logical unit number of data base input file
character*50 cfile
c ! output filename
character*2 coddat(ncode)
c ! allowable letter codes
character*2 code
c ! letter code of data base choice
real delta
c ! angular distance
real dg2rad
c ! pi/180.
character*15 filnam(ncode)
c ! filenames of database
character*40 filnm
c ! output filename
real ggtogc
c ! function
real height
c ! plot height in inches
integer ifor
c ! output format flag
integer iunit
c ! logical unit number of terminal input
real latmax
c ! maximum latitude
real latmin
c ! minimum latitude
real lonmax
c ! maximum longitude
real lonmin
c ! minimum longitude
integer lunit1
c ! logical unit number of line output file
integer lunit2
c ! logical unit number of .mid output file
real rlnmin
c ! minimum longitude
real rlnmax
c ! maximum longitude
real rltmin
c ! minimum latitude
real rltmax
c ! maximum latitude
real rscale
c ! map scale
real scale
c ! plot scale in km-to-inches
real width
c ! plot width in inches
real xmax
c ! greatest x distance (in km) corresponding to lonmax
real xmin
c ! least x distance (in km) corresponding to lonmin
real ymax
c ! greatest y distance (in km) corresponding to latmax
real ymin
c ! least y distance (in km) corresponding to lonmin
parameter (iunit = 5, lunit1 = 21, lunit2 = 22, bunit = 12)
c
data coddat(1), filnam(1) /'c2', 'c2bin.bin'/
data coddat(2), filnam(2) /'co', 'caoutline.bin'/
data coddat(3), filnam(3) /'cf', 'cafaults.bin'/
data coddat(4), filnam(4) /'cl', 'calakes.bin'/
data coddat(5), filnam(5) /'cr', 'caresrvor.bin'/
data coddat(6), filnam(6) /'si', 'scitex.bin'/
data coddat(7), filnam(7) /'pk', 'pkf.bin'/
data coddat(8), filnam(8) /'mf', 'mletal.bin'/
data coddat(9), filnam(9) /'nc', 'norcal.bin'/
data coddat(10), filnam(10) /'bj', 'baja.bin'/
data coddat(11), filnam(11) /'mh', 'mhflts.bin'/
data coddat(12), filnam(12) /'bf', 'banning.bin'/
data coddat(13), filnam(13) /'ap', 'apflts.bin'/
data coddat(14), filnam(14) /'ac', 'apcrp.bin'/
data coddat(15), filnam(15) /'cg', 'coalinga.bin'/
data coddat(16), filnam(16) /'fd', 'socalfold.bin'/
data coddat(17), filnam(17) /'sf', 'sfbayrds.bin'/
data coddat(18), filnam(18) /'nh', 'nevhist.bin'/
data coddat(19), filnam(19) /'nq', 'nevhollp.bin'/
data coddat(20), filnam(20) /'no', 'nevholo.bin'/
data coddat(21), filnam(21) /'np', 'nevlp.bin'/
data coddat(22), filnam(22) /'nd', 'nevundat.bin'/
data coddat(23), filnam(23) /'nu', 'nevundif.bin'/
data coddat(24), filnam(24) /'nv', 'nevvolc.bin'/
data coddat(25), filnam(25) /'st', 'states.bin'/
data coddat(26), filnam(26) /'yl', 'yellowstn.bin'/
data coddat(27), filnam(27) /'so', 'sroutline.bin'/
data coddat(28), filnam(28) /'hw', 'hawaii.bin'/
data coddat(29), filnam(29) /'wm', 'world.bin'/
data coddat(30), filnam(30) /'pt', 'plates.bin'/
data coddat(31), filnam(31) /'cc', 'cacounty.bin'/
data coddat(32), filnam(32) /'f0', 'jen92_hist.bin'/
data coddat(33), filnam(33) /'f1', 'jen92_holo.bin'/
data coddat(34), filnam(34) /'f2', 'jen92_q7.bin'/
data coddat(35), filnam(35) /'f3', 'jen92_q16.bin'/
c
dg2rad = 4.0*atan(1.0)/180.
c
c get latitude, longitude bounds
c
call gtbnds (iunit, ounit, latmin, latmax, lonmin, lonmax)
c
c convert to relative x, y distances in km
c
rltmin = ggtogc(latmin*dg2rad)
rltmax = ggtogc(latmax*dg2rad)
rlnmin = lonmin*dg2rad
rlnmax = lonmax*dg2rad
call refpt ((rltmax + rltmin)/2., (rlnmax + rlnmin)/2.)
call delaz (rltmin, rlnmin, delta, azm, bazm, xmin, ymin)
call delaz (rltmax, rlnmax, delta, azm, bazm, xmax, ymax)
height = 7.
width = 7.
rscale = 62500
c
c set up scaling
c
10 ans = 'd'
call askc (
& 'Enter scaling option: s = map scale, or d = map dimensions',
& ans, ounit)
if (ans .ne. 's' .and. ans .ne. 'd') then
write (ounit, *) 'Please answer "s" or "d": '
goto 10
end if
c
c map scale
c
if (ans .eq. 's') then
20 rscale = askr ('Enter map scale', rscale, ounit)
if (rscale .le. 0.) then
write (ounit, *) 'Scale must be greater than 0; try again'
goto 20
end if
scale = rscale*(2.54e-5)
height = (ymax - ymin)/scale
width = (xmax - xmin)/scale
write (6, 30) height, width
30 format (/, 1x, 'Map will be approximately ', f8.2,
& ' inches high by ', f8.2, ' inches wide')
else
c
c explicit distance scale
c
40 width = askr ('Enter map width (inches)', width, ounit)
if (width .le. 0.) then
write (ounit, *)
& 'Width must be greater than 0; try again'
goto 40
end if
50 height = askr ('Enter map height (inches)', height, ounit)
if (height .le. 0.) then
write (ounit, *)
& 'Height must be greater than 0; try again'
goto 50
end if
scale = max((xmax - xmin)/width, (ymax - ymin)/height)
write (6, 60) nint(scale/2.54e-5)
60 format (/, 'Map will have a scale of 1:', i)
endif
70 ans = 'y'
call askc ('OK (y=yes, n=no, q=quit) ?', ans, ounit)
if (ans .ne. 'y' .and. ans .ne. 'n' .and. ans .ne. 'q') then
write (ounit, *) 'Unrecognized response; try again'
goto 70
end if
if (ans .eq. 'n') then
goto 10
elseif (ans .eq. 'q') then
stop
end if
c
c output format
c
write (ounit, *) ' '
write (ounit, *) 'The following output formats are available'
write (ounit, *) '0 = Default Qplot: 6 coordinate pairs/record,'
1 //' decimal degrees,'
write (ounit, *) ' + west longitude, fortran format (6(f6.4,'
1 //'f7.4)), pen lift=(0.0, 0.0)'
write (ounit, *) '1 = Ascii dump: 1 coordinate pair/record,'
1 //' decimal degrees,'
write (ounit, *) ' - west longitude, fortran format (2(f10.5'
1 //', 5x)), pen lift=(-999.0, -999.0)'
write (ounit, '(a,/)')
1 '2 = MapInfo (.MIF & .MID interchange format)'
75 ifor = 0
ifor = jask ('Enter format [0-2]', ifor)
if (ifor .lt. 0 .or. ifor .gt. 2) then
write (ounit, *) 'Please enter valid format #'
goto 75
end if
c
c open output file
c
80 if (ifor .lt. 2) then
cfile = 'sift.lin'
call askc ('Enter output filename', cfile, ounit)
open (lunit1, file = cfile, status = 'unknown')
else
cfile = 'mapinfo'
call askc ('Enter mapinfo base output filename', cfile, ounit)
do 90 i = 40, 1, -1
if (cfile(i:i) .ne. ' ') then
if (i .gt. 46) then
print *, 'warning: output filenames are being truncated'
i = 46
end if
print *, ' '
filnm = cfile(1:i)//'.mif'
open (lunit1, file = filnm, status = 'unknown')
print *, 'MIF output written to ', filnm
filnm = cfile(1:i)//'.mid'
open (lunit2, file = filnm, status = 'unknown')
print *, 'MID output written to ', filnm
goto 100
end if
90 continue
print *, 'error forming output filname; no blank found'
stop
end if
c
c main database menu
c
100 write (ounit, '(/,a)') 'Data sets available are....'
write (ounit, *) ' ca = California files'
write (ounit, *) ' nv = Nevada files'
write (ounit, *) ' us = Other USA files'
write (ounit, *) ' wd = World files'
write (ounit, *) ' nf = Open new output file'
write (ounit, *) ' ? = This menu'
write (ounit, *) ' x = Exit program'
write (ounit, *) ' '
code = 'ca'
call askc ('Enter code: ', code, ounit)
c
c get the goods
c
if (code .eq. 'x ') then
close (lunit1)
if (ifor .eq. 2) close (lunit2)
stop
else if (code .eq. 'nf') then
close (lunit1)
if (ifor .eq. 2) close (lunit2)
goto 80
else if (code .eq. '?') then
goto 100
else if (code .eq. 'ca') then
200 write (ounit, *) 'California menu'
write (ounit, *) ' co = State outline (1:250,000 scale)'
write (ounit, *) ' f0 = Historic (<200 yrs) faults '//
& '(Jennings, 1992)'
write (ounit, *) ' f1 = Holocene (<10 ka) faults '//
& '(Jennings, 1992)'
write (ounit, *) ' f2 = Quaternary (<700 ka) faults '//
& '(Jennings, 1992)'
write (ounit, *) ' f3 = Quaternary (<1600 ka) faults '//
& '(Jennings, 1992)'
write (ounit, *) ' cl = Major lakes and reservoirs'
write (ounit, *) ' cr = Minor lakes and reservoirs'
write (ounit, *) ' cf = Major Quat. faults (Jennings, 1975)'
c write (ounit, *) ' si = Raster scan of Jennings (1975) Q'//
c & 'uat. flts (slow reading but high detail)'
write (ounit, *) ' cc = County boundaries'
write (ounit, *) ' mf = Mammoth Lakes roads, geog'//
& 'raphy'
c write (ounit, *) ' mh = Morgan Hill region faults (Herd)'
c write (ounit, *) ' bf = Banning region faults (Matti)'
c write (ounit, *) ' nc = N. CA faults (see Castillo/Ellswr'//
c write (ounit, *) ' ap = Alquist-Priola flts in SF Bay reg'//
c & 'ion'
write (ounit, *) ' ac = Active creep sites in SF Bay region'
write (ounit, *) ' cg = Coalinga anticline'
write (ounit, *) ' fd = S. Ca. Cenoz. anticlinal fold axes'
write (ounit, *) ' sf = Major highways in SF Bay region'
write (ounit, *) ' ? = This menu'
write (ounit, *) ' x = Exit California menu'
write (ounit, *) ' '
code = 'co'
210 call askc ('Enter code: ', code, ounit)
if (code .eq. 'x') then
goto 100
elseif (code .eq. '?') then
goto 200
else
call rdcode (code, icode, coddat, ncode)
if (icode .eq. 0) then
print *, 'Unrecognized code; try again'
goto 210
else
c
c exception to handle copyrighted CDMG fault files
c
if (icode .gt. 31 .and. icode .lt. 36) then
call rdfile (lunit1, latmax, latmin, lonmax, lonmin,
& scale, prefix2//filnam(icode), bunit, ifor, lunit2)
else
call rdfile (lunit1, latmax, latmin, lonmax, lonmin,
& scale, prefix1//filnam(icode), bunit, ifor, lunit2)
endif
end if
end if
goto 210
c
c Nevada files
c
else if (code .eq. 'nv') then
300 write (ounit, *) 'Nevada menu (* => recommended)'
write (ounit, *) ' nh = historical faults (*)'
write (ounit, *) ' nq = Holocene-Late Pleist. faults (*)'
write (ounit, *) ' no = Holocene faults (*)'
write (ounit, *) ' np = Late Pleist. faults (*)'
write (ounit, *) ' nd = undated faults (huge)'
write (ounit, *) ' nu = undifferentiated faults'
write (ounit, *) ' nv = faults offsetting volcanics'
write (ounit, *) ' ? = This menu'
write (ounit, *) ' x = Exit Nevada menu'
write (ounit, *) ' '
code = 'nh'
310 call askc ('Enter code: ', code, ounit)
if (code .eq. 'x') then
goto 100
elseif (code .eq. '?') then
goto 300
else
call rdcode (code, icode, coddat, ncode)
if (icode .eq. 0) then
print *, 'Unrecognized code; try again'
goto 310
else
call rdfile (lunit1, latmax, latmin, lonmax, lonmin,
& scale, prefix1//filnam(icode), bunit, ifor, lunit2)
end if
end if
goto 310
c
c Rest of USA
c
else if (code .eq. 'us') then
400 write (ounit, *) 'USA menu'
write (ounit, *) ' st = U.S. state boundaries'
write (ounit, *) ' yl = Yellowstone faults, boundaries,'//
& ' lakes, calderas'
write (ounit, *) ' so = Snake River Plain outline'
write (ounit, *) ' hw = Hawaii elev, vents, calderas, etc.'
write (ounit, *) ' ? = This menu'
write (ounit, *) ' x = Exit USA menu'
write (ounit, *) ' '
code = 'st'
410 call askc ('Enter code: ', code, ounit)
if (code .eq. 'x') then
goto 100
elseif (code .eq. '?') then
goto 400
else
call rdcode (code, icode, coddat, ncode)
if (icode .eq. 0) then
print *, 'Unrecognized code; try again'
goto 410
else
call rdfile (lunit1, latmax, latmin, lonmax, lonmin,
& scale, prefix1//filnam(icode), bunit, ifor, lunit2)
end if
end if
goto 410
c
c World
c
else if (code .eq. 'wd') then
500 write (ounit, *) 'World menu'
write (ounit, *) ' wm = World map'
write (ounit, *) ' pt = World tectonic plate boundaries'
write (ounit, *) ' bj = Baja Mexico faults'
write (ounit, *) ' ? = This menu'
write (ounit, *) ' x = Exit World menu'
write (ounit, *) ' '
code = 'wm'
510 call askc ('Enter code: ', code, ounit)
if (code .eq. 'x') then
goto 100
elseif (code .eq. '?') then
goto 500
else
call rdcode (code, icode, coddat, ncode)
if (icode .eq. 0) then
print *, 'Unrecognized code; try again'
goto 510
else
call rdfile (lunit1, latmax, latmin, lonmax, lonmin,
& scale, prefix1//filnam(icode), bunit, ifor, lunit2)
end if
end if
goto 510
c
c test module
c
else if (code .eq. 'c2') then
call rdfile (lunit1, latmax, latmin, lonmax, lonmin,
& scale, filnam(1), bunit, ifor, lunit2)
goto 100
c
c last test of main menu
c
else
write (ounit, *) 'Unrecognized code; try again: '
goto 100
endif
c
end