program loadtopo c c Program by John Boursy, November 1982. c c This program is designed to read the topographic data base in c its 4812-character per record format, compress it into blocks c of 30 minutes by 30 minutes, and write out the compressed c records. c c The compressed records are written out to separate files, c according to the longitude degrees. The longitudes handled c range from 66 to 125 degrees West Longitude. The compressed c records are written to files which have logical unit numbers of c the longitude degrees minus 46. c c Another program, mergetopo, reads all of these intermediate files c and merges them into a single output file. c c ****************************************************************** c c Following is a description of the input topographic data: c c The topographic data was obtained from the National Geophysical c and Solar-Terrestrial Data Center in Boulder, Colorado. The c mailing address is: c c U.S. Dept. of Commerce c NOAA/EDIS/NGSDC (D623) c 325 Broadway c Boulder, Colorado 80303 c c The topographic information was obtained on four magnetic tapes. c They are 9-track, 1600 bpi, EBCDIC, without labels, and without c serial numbers. There is one record per block, with each record c 4812 characters long. The tapes covered the following bands of c longitude, and contained the following number of records: c c Tape No. Longitude Band No. of Records c c 1 66-86 deg W 3960 c 2 86-100 deg W 4092 c 3 100-109 deg W 2340 c 4 109-125 deg W 3876 c c The information which came with the tapes indicated the following c people to contact: c c (a) Format/tape questions -- David Dater (303) 497-6120 c c (b) Other questions -- David Clark (303) 497-6123 (Actually, it c was Allen Hittleman, but he has since left NGDSC, and we talk c with Dave Clark.) c c To run this program, LOADTOPO, we copied the data from the four c tapes into four disk files, named TOPO1.DAT, TOPO2.DAT, TOPO3.DAT, c and TOPO4.DAT. During this copying, we converted from EBCDIC to c ASCII; this program does NOT do character transliteration. c c Once we have the four disk files, we run LOADTOPO four times, c once for each disk file. The name of the input file is entered c on a data card that is read by LOADTOPO. LOADTOPO handles the c proper generation of the output files by using the longitude of c the data; it does not depend on the name of the input files. c c ****************************************************************** c c Following is a discussion of the format of the 4812-character c records in the input files: c c The data on the tape is 30-second point elevation data. Each c 4812-character record contains a 5-minute latitude by one degree c of longitude series of elevation records. The first 12 characters c define the southeast corner of this geographical region and the c remaining 4800 characters represent a sequence of 1200 elevation c values, each being four characters long. c c The sequencing of records is along latitude bands of 30-seconds c which index from 0'0" to 59'30" in longitude spacing. Thus, there c are 120 elevation values per 30-second latitude band. c c *************** ... *************** 5'00" c * * * * * * c * 1200 * 1199 * * 1082 * 1081 * c * * * * * * c *************** ... *************** c . . . . c . . . . c . . . . c *************** ... *************** c * * * * * * c * 240 * 239 * * 122 * 121 * c * * * * * * c *************** ... *************** c * * * * * * c * 120 * 119 * * 2 * 1 * c * * * * * * c *************** ... *************** 0'00" c 60'00" 1'00" 0'00" c c The format may be described as follows: c c Character FORTRAN Description c Position FORMAT c c 1-4 I4 Latitude degrees of southeast c corner of geographical area. c (Note: characters 1 and 2 are c blank.) c c 5-8 I4 Latitude minutes of southeast c corner of geographical area. c (Note: characters 5 and 6 are c blank.) c c 9-12 I4 Longitude degrees of southeast c corner of geographical area. c (Note: characters 9 and 10 c may be blank. All values are c unsigned west longitude.) c c 13-4812 1200 I4 Elevation in meters above mean c sea level of the 1200 c elevations. c c Each sequence of 12 records represents a 1 degree by 1 degree c area. c c ****************************************************************** c c The compressed records will be variable length. The more c compression, the shorter they will be. The records will be c composed of an array of integers. (Integer*4 integers will be c used throughout.) Here's a description of the record: c c 1st word -- This is the key. It is composed of longitude and c latitude of the southeast corner of the 30 minute c by 30 minute block. North latitude is positive, c and south latitude is negative. West longitude is c positive, and east longitude is negative. The c longitude will be multiplied by 100,000 and added to c ten times the latitude. For example, if the c southeast corner is 32.5 degrees North latitude, c and 109 degrees West longitude, the key would be c an integer with the value of 10900325. c c 2nd word -- This contains various descriptive information about c this record. This is broken up into groups of c bits. Since this is for a VAX 11/780, each integer c contains 32 bits. In keeping with VAX conventions, c the least significant bit is bit 0, while the most c significant bit is bit 31. (This is just the c opposite of Univac 1100 and Honeywell 6000, so c watch out in dealing with this.) Here is how the c bits are used: c c bit range description c c 0-3 contains the number of bits c per elevation for each c elevation in this record. c 4-10 contains the value by which c the compressed value must be c multiplied to get the actual c value. Initially, this will c be 10 because all elevations c are divisible by 10; we c divide them by 10, store c them, and then multiply by 10 c and then retrieve them. c 11-15 Unused. c 16 Indicates whether the values c in this record are in meters c or feet. This is 1 if the c values are in meters; 0 if c values are in feet. c 17 Indicates whether some c elevations in this record are c missing. This is 1 if some c are missing; otherwise, 0. c 18-31 Unused. c c 3rd word -- contains the bias. This value must be added to all c values that are stored to recreate the original c values. The stored values must be multiplied by c the value in bits 4 through 10 of the second word; c then, the bias is added. c c 4th word -- contains the number of words after this word. The c words after this contain the compressed elevations. c If this is zero, there are not words after this, c and this means that the entire 30 minute by 30 c minute block is in the ocean with a height above c mean sea level of zero. c c 5th word c thru end -- contains the compressed elevations. Each word c contains several elevations. The number per word c is equal to 32 divided by the number of bits per c elevation given in bits 0 thru 4 of the second c word. If this division ends up with a remainder, c the extra bits are simply set to zero. We do not c try to store elevations across word boundaries. c For example, if there are 10 bits per elevation c (according to bits 0 thru 4 of the second word), c then bits 0 thru 9 are used for the first c elevation, bits 10 thru 19 are used for the second c elevation, bits 20 thru 29 are used for the third c elevation, and bits 30 and 31 are set to zero. c c ****************************************************************** c c Next is a description of some adjustments which we make to the c topographic data base before we write out our files. These are to c take account of some problems with the source data. They are: c c 1. Although the entire topographic data base is supposed to be in c meters, there are a couple of areas in which the values are in c feet. As noted above, we can use bit 16 of the second word of c the output record to keep track of whether we are using meters c or feet. We must set this bit to zero when we are using feet, c in the few areas which have feet rather than meters. These c areas are: c c a. The 1-degree by 1-degree block between 44 and 45 degrees c North latitude and 87 and 88 degrees West longitude. c c b. The 2-degree by 1-degree block between 25 and 27 degrees c North latitude and 97 and 98 degrees West longitude. c c 2. There are also some areas where the elevation for a lake is c given as zero meters, when it should be 180 meters. We adjust c the zero elevations to 180 before writing out the records. c These areas are: c c a. The 2-degree by 30-minute block between 42 and 44 degrees c North latitude and between 87 degrees and 87 degrees 30 c minutes West longitude. c c b. The 1-degree by 30-minute block between 42 degrees 30 c minutes and 43 degrees 30 minutes North latitude and c between 86 degrees 30 minutes and 87 degrees West c longitude. c c c. The 30-minute by 30-minute block between 44 degrees 30 c minutes and 45 degrees North latitude and between 86 c degrees 30 minutes and 87 degrees West longitude. c c ****************************************************************** c c The following statement is the first statement. c c ****************************************************************** c parameter (maxwords=1800) ! max # of output words dimension ielev(1200) dimension jelev(120,120) dimension latd(12),latm(12),lond(12) character*4 celev(1200) data lutopo/11/ integer out integer outtemp/12/ ! temporary output file dimension lulist(60) ! intermediate output files data lulist/60*0/ data numlu/0/ ! # of intermediate output files integer pout/6/ data lonold/-300/ dimension iout(maxwords) dimension ibitcnt(32) data ibitcnt/32*0/ logical first79 /.true./ character*30 topofile data in/5/ c c ****************************************************************** c c Following is the statement function section. c c ****************************************************************** c alog2(x)=log(x)/log(2.) c This computes the logarithm of x to the base 2. c c ****************************************************************** c c The following statement is the first executable statement. c c ****************************************************************** c call ptime (0,pout) c write (pout,850) 850 format (/,'$File name of uncompressed file: ') read (in,851) topofile 851 format (a30) if (topofile.eq.' ') go to 9000 write (pout,852) topofile 852 format ('0The uncompressed input file is ',a30) c write (pout,814) 814 format ('0Welcome to Loadtopo',//,10x,'Key',8x,'Desc',2x,'Bias', 2 3x,'Bits',1x,'# wds',3x,'Max',3x,'Min') c open (unit=lutopo,status='old',file=topofile, 2 form='formatted',readonly,shared,iostat=iostat,err=2000) open (unit=outtemp,status='scratch',form='unformatted',recl=1807, 2 recordtype='variable',iostat=iostat,err=3100) c i=0 ! counts input records j=0 ! counts output records ifunny=0 ! counts # of missing elevations none=0 ! counts # of output records with missing elevations jdelta=10 ! resolution of points minall=10000000 maxall=-minall minlat=100 maxlat=-minlat minlon=200 maxlon=-minlon minwd=maxwords*2 ! size of minimum compressed output record maxwd=-minwd ! size of maximum compressed output record mindel=10000000 maxdel=-mindel rewind lutopo itemp=0 ! keeps track of # of records on temporary file c 10 continue do 100 loop=1,12,1 read (lutopo,801,end=1000) latd(loop),latm(loop),lond(loop), 2 (celev(loop2),loop2=1,1200,1) 801 format (3i4,1200a4) i=i+1 if (latd(loop).lt.minlat) minlat=latd(loop) if (latd(loop).gt.maxlat) maxlat=latd(loop) if (lond(loop).lt.minlon) minlon=lond(loop) if (lond(loop).gt.maxlon) maxlon=lond(loop) c do 20 loop2=1,1200,1 if (celev(loop2).eq.' ') then ielev(loop2)=-9999 ! set to find missing elev. later else read (celev(loop2),802) ielev(loop2) 802 format (i4) if (mod(ielev(loop2),10).ne.0) write (pout,803) 2 latd(loop),latm(loop),lond(loop),ielev(loop2) 803 format ('0Not divisible by 10 at',3i6,'; elevation is',i6) endif 20 continue c do 25 loop2=1,10,1 istart=(loop2-1)*120+1 iend=istart+119 iy=(loop-1)*10+loop2 do 23 loop3=istart,iend,1 jelev(loop3-istart+1,iy)=ielev(loop3) 23 continue 25 continue c 100 continue c c ****************************************************************** c c In this next section, we make the adjustments necessary to c compensate for incorrect data in the input data base: c c ****************************************************************** c if ((latd(1).eq.42.or.latd(1).eq.43).and.lond(1).eq.87) then c We adjust in the range from 42-00 to 44-00 and 87-00 to 87-30 c by changing zero elevations to elevations of 180 meters. c do 60 loop=1,60,1 do 50 loop2=1,120,1 if (jelev(loop,loop2).eq.0) jelev(loop,loop2)=180 50 continue 60 continue c endif c if ((latd(1).eq.42.or.latd(1).eq.44).and.lond(1).eq.86) then c We adjust in the range from 42-30 to 43 and 86-30 to 87-00 or c in the range from 44-40 to 45-00 and 86-30 to 87-00 by changing c zero elevations to elevations of 180 meters. c do 80 loop=61,120,1 do 70 loop2=61,120,1 if (jelev(loop,loop2).eq.0) jelev(loop,loop2)=180 70 continue 80 continue c endif c if (latd(1).eq.43.and.lond(1).eq.86) then c We adjust in the range from 43-00 to 43-30 and 86-30 to 87-00 c by changing the zero elevations to elevations of 180 meters. c do 95 loop=61,120,1 do 90 loop2=1,60,1 if (jelev(loop,loop2).eq.0) jelev(loop,loop2)=180 90 continue 95 continue c endif c min=10000000 max=-min c do 120 loop=1,120,1 do 110 loop2=1,120,1 if (jelev(loop,loop2).lt.min) min=jelev(loop,loop2) if (jelev(loop,loop2).gt.max) max=jelev(loop,loop2) 110 continue 120 continue c if (min.lt.minall) minall=min if (max.gt.maxall) maxall=max c c ****************************************************************** c c Next, we determine if we should have our elevations in meters or c in feet. c c ****************************************************************** c imetric=1 ! initialize to show elevations in meters if ((latd(1).eq.25.or.latd(1).eq.26).and.lond(1).eq.97) imetric=0 if (latd(1).eq.44.and.lond(1).eq.87) imetric=0 c if (lond(1).eq.79.and.first79) then c c This next section does some preliminary loading for some 30 c minute by 30 minute blocks that are missing from the sequential c topo data base file. These are all in the ocean (so they will c be loaded with zeros), but they are close enough that they are c considered in doing 2 to 10 mile profiles. c first79=.false. ! go thru here no more than once max=0 min=0 iword=0 iout(1)=0 ibits=0 idesc=0 call mvbits (ibits,0,4,idesc,0) call mvbits (jdelta,0,7,idesc,4) call mvbits (imetric,0,1,idesc,16) ibias=-jdelta numlu=numlu+1 lulist(numlu)=79-46 c out=lulist(numlu) open (unit=out,status='new',form='unformatted',recl=1805, 2 recordtype='variable',iostat=iostat,err=3000) c c ****************************************************************** c c In this next section, we create some compressed records with c elevations of 0 above mean sea level. We have missing stuff on c the input file in these areas of the sea near land. c c ****************************************************************** c do 5 loop=250,265,5 ! from 25 deg to 26 degrees, 30 minutes none=none+1 key=7950000+loop ! only at 79 degrees, 30 minutes write (79-46,iostat=iostat,err=4100) key,idesc,ibias,iword, 2 (iout(loop2),loop2=1,iword,1) write (pout,804) key,idesc,ibias,ibits,iword,max,min j=j+1 5 continue c endif c do 30 loop=1,numlu,1 if (lulist(loop)+46.eq.lond(1)) go to 40 30 continue loop=0 c 40 continue if (loop.eq.0) then numlu=numlu+1 lulist(numlu)=lond(1)-46 out=lulist(numlu) open (unit=out,status='new',form='unformatted',recl=1805, 2 recordtype='variable',iostat=iostat,err=3000) else out=lulist(loop) endif c if (lond(1).ne.lonold.and.lonold.ne.-300) then itemp=0 ! counts records being moved from temp file rewind outtemp 600 continue read (outtemp,end=700) key,max,min,idesc,ibias,iword, 2 (iout(loop5),loop5=1,iword,1) write (lonold-46,iostat=iostat,err=4000) key,idesc,ibias,iword, 2 (iout(loop5),loop5=1,iword,1) ibits=0 call mvbits (idesc,0,4,ibits,0) write (pout,804) key,idesc,ibias,ibits,iword,max,min 804 format (1x,i14,2x,z8,5i6) go to 600 700 continue rewind outtemp endif lonold=lond(1) c do 125 loop=1,12,1 if (latd(1).ne.latd(loop)) write (pout,815) latd(1),loop, 2 latd(loop),loop,lond(loop) 815 format ('0*** Non-matching Latitude: Latd(1) is',i4, 2 '; Latd(',i2,') is',i4,'; Lond(',i2,') is',i4) if (lond(1).ne.lond(loop))write (pout,816) lond(1),loop, 2 lond(loop),loop,latd(loop) 816 format ('0*** Non-matching Longitude: Lond(1) is',i4, 2 '; Lond(',i2,') is',i4,'; Latd(',i2,') is',i4) if (latm(loop).ne.(loop-1)*5) write (pout,817) loop,latm(loop), 2 loop,latd(loop),loop,lond(loop) 817 format ('0*** Bad Latitude Minutes: Latm(',i2,') is',i4, 2 '; Latd(',i2,') is',i4,'; Lond(',i2,') is',i4) 125 continue c do 500 loop=1,2,1 do 400 loop2=1,2,1 min=10000000 max=-min istartx=1+(loop-1)*60 istarty=1+(loop2-1)*60 iendx=istartx+59 iendy=istarty+59 izero=0 c do 200 loop3=istartx,iendx,1 do 150 loop4=istarty,iendy,1 if (jelev(loop3,loop4).lt.min.and.jelev(loop3,loop4).ne.-9999) 2 min=jelev(loop3,loop4) if (jelev(loop3,loop4).gt.max.and.jelev(loop3,loop4).ne.-9999) 2 max=jelev(loop3,loop4) if (jelev(loop3,loop4).eq.0) izero=izero+1 150 continue 200 continue c lat=latd(1)*10 lon=lond(1)*10 if (loop.eq.2) then lon=lon+5 luout=outtemp itemp=itemp+1 else luout=out endif if (loop2.eq.2) lat=lat+5 key=lon*10000+lat idesc=0 idelta=max-min if (idelta.lt.mindel) mindel=idelta if (idelta.gt.maxdel) maxdel=idelta c ibias=min-jdelta c By setting the bias to be min-jdelta, we will always have c positive values of compressed elevations. This allows us to c use zero for the undefined values, and prevents having to deal c with negative values at all. The lowest possible value that can c be stored, then, is min-jdelta. But since the actual lowest c value is min, we are using the value of min-jdelta to represent c the situation where the elevation is missing altogether. c idelta=(idelta+jdelta)/jdelta c We set up the new idelta to be, basically, the old idelta divided c by jdelta. However, we actually use idelta+jdelta in the c numerator because we have to allow for the fact that we are c really beginning at min-jdelta instead of min. The new idelta is c used to determine the range of values which we have to store, c taking into account that they will be multiplied by jdelta on c retrieval from the compressed data base. c if (izero.eq.3600.or.idelta.eq.0) then if (izero.ne.3600) write (pout,819) max,min 819 format (5x,'All same elevation; min =',i5,'; max =',i5) iword=0 none=none+1 iout(1)=0 ibits=0 call mvbits (ibits,0,4,idesc,0) ! # of bits/elevation call mvbits (jdelta,0,7,idesc,4) ! set bias call mvbits (imetric,0,1,idesc,16) ! set it for metric or not c if (luout.eq.outtemp) then ! write to temp file write (luout,iostat=iostat,err=4000) key,max,min,idesc, 2 ibias,iword,(iout(loop5),loop5=1,iword,1) else write (luout,iostat=iostat,err=4100) key,idesc,ibias,iword, 2 (iout(loop5),loop5=1,iword,1) write (pout,804) key,idesc,ibias,ibits,iword,max,min endif c j=j+1 else ibits=alog2(float(idelta)) if (idelta.ge.2**ibits) ibits=ibits+1 if (ibits.lt.1) ibits=1 ibitcnt(ibits)=ibitcnt(ibits)+1 iperwd=32/ibits ! # of compressed elevations/word ipoint=0 ! where in word it is going; ipoint=2 is 2nd in wd iword=1 jfunny=0 ! # of missing elevations in this block maxval=2**ibits ! maximum possible value to be stored c do 300 loop3=istartx,iendx,1 do 250 loop4=istarty,iendy,1 if (jelev(loop3,loop4).eq.-9999) then jelev(loop3,loop4)=0 jfunny=jfunny+1 else jelev(loop3,loop4)=(jelev(loop3,loop4)-ibias)/jdelta if (jelev(loop3,loop4).gt.maxval) write (pout,821) 2 jelev(loop3,loop4),maxval 821 format (' *** Truncation in writing value of',i6, 2 '; maximum allowed value is',i6) endif ipoint=ipoint+1 if (ipoint.gt.iperwd) then ipoint=1 iword=iword+1 iout(iword)=0 if (iword.gt.maxwords) write (pout,813) key 813 format ('0*** Too many output words for key = ',i12) endif c ifirst=0 call mvbits(jelev(loop3,loop4),ifirst,ibits,iout(iword), 2 (ipoint-1)*ibits) 250 continue 300 continue ifunny=ifunny+jfunny c call mvbits (ibits,0,4,idesc,0) ! set # bits/elevation call mvbits (jdelta,0,7,idesc,4) ! set bias call mvbits (imetric,0,1,idesc,16) ! set it as metric or not if (jfunny.gt.0) then missing=1 else missing=0 endif call mvbits (missing,0,1,idesc,17) ! set whether missing elev c if (iword.lt.minwd) minwd=iword if (iword.gt.maxwd) maxwd=iword c if (luout.eq.outtemp) then ! write to temp file write (luout,iostat=iostat,err=4000) key,max,min,idesc, 2 ibias,iword,(iout(loop5),loop5=1,iword,1) else write (luout,iostat=iostat,err=4100) key,idesc,ibias,iword, 2 (iout(loop5),loop5=1,iword,1) write (pout,804) key,idesc,ibias,ibits,iword,max,min endif c j=j+1 endif 400 continue 500 continue c go to 10 c 1000 continue c We come through here when we hit the end-of-file on input. c if (itemp.gt.0) then ! move from temp file to output file rewind outtemp 1200 continue read (outtemp,end=1300) key,max,min,idesc,ibias,iword, 2 (iout(loop5),loop5=1,iword,1) write (out,iostat=iostat,err=4000) key,idesc,ibias,iword, 2 (iout(loop5),loop5=1,iword,1) ibits=0 call mvbits (idesc,0,4,ibits,0) write (pout,804) key,idesc,ibias,ibits,iword,max,min go to 1200 1300 continue endif c write (pout,808) i,j 808 format (i12,' input records;',i12,' output records') write (pout,809) none 809 format (i12,' records that are all zeros') write (pout,812) ifunny 812 format (i12,' missing elevations') write (pout,805) minall,maxall 805 format ('0Overall minimum elevation='i5,/, 2 '0Overall maximum elevation=',i5) write (pout,806) minlat,maxlat,minlon,maxlon 806 format ('0Overall minimum latitude degrees is',i6,/, 2 '0Overall maximum latitude degrees is',i6,/, 3 '0Overall minimum longitude degrees is',i6,/, 4 '0Overall maximum longitude degrees is',i6) write (pout,818) minwd,maxwd 818 format ('0Overall minimum # of compressed output words =',i5,/, 2 '0Overall maximum # of compressed output words =',i5) c do 950 loop=1,32,1 if (ibitcnt(loop).ne.0) write (pout,820) loop,ibitcnt(loop) 820 format (' With',i3,' bits/elevation:',i6,' Records') 950 continue c go to 9000 ! closing up everything c 2000 continue c We come through here if we have trouble opening the input file. write (pout,807) iostat 807 format ('0*** Trouble in opening TOPO1.DAT; iostat is',i4,' ***') go to 9000 ! closing up everything c 3000 continue c We come through here if we have trouble opening an intermediate c output file. write (pout,810) out,iostat 810 format ('0*** Trouble in opening logical unit',i3, 2 '; iostat is',i4,' ***') go to 9000 ! closing up everything c 3100 continue c We come through here if we have trouble opening our temporary c file, outtemp. write (pout,810) outtemp,iostat go to 9000 ! closing up everything c 4000 continue c We come through here if we have an error in doing our c non-formatted write. write (pout,811) out,iostat,key,idesc,ibias,iword 811 format ('0*** Error in writing to logical unit',i3, 2 '; iostat is',i4,' ***',/,'0Key is',i12,'; Idesc is ',z8, 3 '; Ibias is',i6,'; Iword is',i6) go to 9000 ! closing up everything c 4100 continue c We come through here if we have an error in doing our c non-formatted write. write (pout,811) luout,iostat,key,idesc,ibias,iword c 9000 continue c We come through here to close up everything call ptime(1,pout) stop end