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