The code excerpt below converts data from widely-used ARC-MOSS format into CAMRIS format. It included coordinate transformations to and from lat.-long. and the Lambert conical projection used in geographic survey maps. Some elementary calculus is used to invert the power series expressing the transformation.
c MOSVEC.FOR 11/1/94 DLK c moss to CAMRIS format conversion, including c conversion routines between Lambert (X,Y meters) coords and lat.-lon. c for Costa Rica Norte (Ocotepeque) Chart (Clark Spheroid 1866). Real*8 x,y, lat, lon Character*40 afil 1 print *,'Enter file name' read(*,101) afil 101 format(a40) if(afil.eq.' ') stop open(unit=1, file=afil) open(unit=2, file=afil(1:5)//'.tst') 2 read(1,*, end=3) id, idum,nvert if(nvert.le.18) write(2,102) nvert,id,idum write(*,102) nvert,id,idum 102 format(i5,15x,i5,5x,i5) do i=1,nvert read(1,*) x,y call lamltln(x,y,lat,lon) if(nvert.le.18) write(2,103) lon,lat c write(*,104) id,nvert,x,y,lat,lon 103 format(6f13.6) c104 format(2i5,4f13.6) enddo goto 2 3 close(1) close(2) stop END c********************************* Subroutine Lamltln(x,y, lat,lon) Real*8 x,y, lat, lon c x,y: lambert coordinates in meters; lat, lon: in degrees. c converts x,y, to lat,lon unless x=0, when lat,lon converted to x,y. real*8 dyp common/newton/ dyp c dyp: d(Yprime)/d(dltlat), from Yprime function, for Newton's method c convergence on inverse of Yprime(dltlat). c Misc. variables real*8 dltlat, dltlon, lat0, lon0, xp, yp, ypp, theta, R real*8 epsilon, ypnew, rpd Real*8 Yprime c Chart parameters for Costa Rica Norte Real*8 lat0d, lat0m, lon0d, lon0m, fe, fn, r0fn Real*8 minusa(6) data lat0d,lat0m, lon0d,lon0m/ 10,28, 84,20/ data fe, fn, r0fn/ 5.D5, 2.71820522D5, 3.48D7/ data minusa/ +1843.3251238D0,.000972746D0,2.631931D-5,3.497D-10,5.64D-13,2.D-17/ rpd=3.141592653589793D0/180.D0 c radians per degree, for trig routines... lat0=lat0d*60. + lat0m lon0=lon0d*60. + lon0m c lat0, lon0: in minutes. if(x.ne.0) go to 1 c Conversion from lat-lon to lambert. dltlat=lat*60. - lat0 dltlon= - (lon*60. - lon0) c so that positive dltlon means EAST of lon0...!...... yp=Yprime(dltlat) R=r0fn-yp theta=dltlon*dsin(rpd*lat0/60.) xp=R*dsin(rpd*theta/60.) ypp=xp*dtan(rpd*theta/120.) x=fe+xp y=yp+ypp c print *, lat, lon, lat0, lon0, dltlat, dltlon c print *, yp, R, theta c print *, xp, ypp, x, y Return c Conversion from Lambert to Lat.-Lon. c Yprime(dltlat) is inverted using convergence algorithm. 1 xp=x-fe theta=datan(xp/(r0fn-y))*60./rpd ypp=xp*dtan(rpd*theta/120.) yp=y-ypp c Converge on dltlat from y prime. Invert linear term as initial guess. dltlat=(yp-fn)/minusa(1) niter=0 2 ypnew=Yprime(dltlat) epsilon=ypnew-yp niter=niter+1 dltlat=dltlat - epsilon/dyp if(epsilon.eq.0) go to 4 if(Dabs(epsilon).lt.1.D-12*yp .and. niter.gt.4) go to 4 if(niter.lt.20) go to 2 4 dltlon=theta/dsin(rpd*lat0/60.) c print *,niter, epsilon, yp, dltlat lat=(dltlat+lat0)/60. lon=( - dltlon+lon0)/60. c (Again, positive dltlon means East of lon0...) END c*********************** Real*8 Function Yprime(dltlat) real*8 dltlat real*8 dyp common/newton/ dyp c Side effect: computes d(Yprime)/d(dltlat), for Newton's method c convergence on inverse function of Yprime. Real*8 dlp c Chart parameters for Costa Rica Norte Real*8 fn Real*8 minusa(6) data fn/2.71820522D5/ data minusa/ +1843.3251238D0,.000972746D0,2.631931D-5,3.497D-10,5.64D-13,2.D-17/ Yprime=fn dyp=minusa(1) dlp=1. do i=1,6 dlp=dlp*dltlat Yprime=Yprime + minusa(i)*dlp if(i.ne.6) dyp=dyp + (i+1)*minusa(i+1)*dlp enddo return END |
The following code was used to create WordPerfect Graphic format files from CAMRIS graphic data, and included run-length data compression.
I also implemented linear and harmonic stretching of map figures and
overlays, for reconciling measurement and digitization errors (not shown).
c CAMWPG.FOR 10/20/95 dlk c CAMRIS--.WPG image conversion routine (excerpts) c this version sets color palette to CAMRIS's, ignores bright/cont, c For color wpg's (even if wp5.0 gets its own spec wrong, other progs don't). c. c.... [some initialization code skipped] c. c****Begin writing .WPG output file. c .WPG file prefix (16 bytes--never varies) call wb (255) a40='WPC' call wst (a40) i4=16 call wl (i4) c NB error in spec: "Start of Data" is long, not short call wb (1) call wb (22) call wb (1) call wb (0) call ws (0) call ws (0) c Data Records call wb (15) call wb (6) c "Start wpg data" record call wb (1) call wb (0) call ws ((1200/ires)*nx) call ws ((1200/ires)*ny) c Image size in wpu (WordPerf units: 1/1200ths of an inch) c (Numbers should be such that there is no roundoff error c here--printer resolutions are usually sub-multiples of c 1200, so there should be no problem). C Set wpg's color map to CAMRIS palette. call wb(14) call wb(255) call ws(4+3*256) c write Color map record call ws(0) call ws(256) do i=0,255 do j=1,3 call wb(ipal(i,j)) enddo enddo call wb (11) c write Bitmap type 1 record c Following 5 bytes before preamble are for record length in longint form call wb (255) lenpos=ipos c save position of length longint for later insertion i4=0 call wl(i4) c fill longint with 0 for now c bitmap preamble--10 bytes call ws (nx) call ws (ny) c nx X ny pixels call ws (idepth) c color depth, in bits (use 1 or 8) call ws (ires) call ws (ires) c set "bitmap source resolution" to desired printed resolution c (see CAMWPG.DAT) c read/write bitmap data call initq Do jy=0,ny-1 Do jx=0,nx-1 Call Pixval(jx,jy,icolor) c Call inqclr(jx,jy,icolor) c (faster--why not?) Call Setcol(256-icolor) c print *,'c x y:',icolor,jx,jy Call Ptabs(jx,jy) c 'invert' screen image colors to show .WPG write progress. if(idepth.eq.8) then C PL version 5/19 use CAMRIS clr index directly--skip other stuff below... iclr=icolor Cc 8-bit color index translation... Cc For some reason our system (wp 5.1, hp 550c) always prints using a Cc uniform greyscale palette, no matter what you tell it, 0=blk, 255=wht... Cc just averaging what I guess are RGB intensities here... C Cc iclr=(ipal(icolor,1)+ipal(icolor,2)+ipal(icolor,3))/3. + .5 Cc iclr=ipal(icolor,1)+ipal(icolor,2)+ipal(icolor,3) C Cc Bright/contrast adjustments. bright=0(black) to 1 Cc contr=0 for no contrast, 1 normal, >1 high (will saturate). Cc set rdp+grp+blp=1, rdp,blp,grp>=0 -- how much r,g,b, resp., contribute to Cc brightness (i.e., set grp>blp>rdp) C C clr=4.*( ipal(icolor,1)*rdp + ipal(icolor,2)*grp + C + ipal(icolor,3)*blp ) - 127.5 C clr=clr*contr + (bright*255) C iclr=max(0,min(255, clr+.5 )) call qbyte(iclr) else ibitpos=7-Mod(jx,8) c Highbit of byte first (dumb). if(ibitpos.eq.7) ibyte=0 c Set color bit (this is the place to tweak 1-bit translation, if needed) if(icolor.eq.22) ibyte=ibyte+2**ibitpos c (0=black, 22=white, in CAMRIS color map) c 0=black, 1=white in .WPG bitmap if(ibitpos.eq.0 .or. jx.eq.nx-1) call qbyte(ibyte) endif enddo call flushq enddo len=ipos-(lenpos+4) c Save bitmap length to insert in start of record call wb (16) call wb (0) c end WPG data record c close output file, and reopen in random access, to stuff in bitmap length close(1) open(unit=1, file=awpg, form='binary', access='direct', + blocksize=512,recl=1) i4 = len c length of record (longint format) i3=i3+128 c Set highest bit of long integer for longint signal write(1,rec=lenpos) i2 c The high word goes first into file (so signal shows). write(1,rec=lenpos+2) i0 c followed by low word. close(1) write(*,*) nx,ny Pause 'Done. Press Enter...' Call Closeg Stop ' ' END c***********end Main*************************************** c*********Queue bytes for compressed output to file******** Subroutine qbyte(ibyte) integer iqueue(0:126),nrep,nqued c print *,'q: start nr,nq:',nrep,nqued c pause ' ' if(nqued.ne.0) go to 2 1 iqueue(0)=ibyte nqued=1 nrep=1 c print *,'q: 1st n,nq:',nrep,nqued c pause ' ' return 2 if(nrep.lt.3) then c in non-repeat state (less than 3 byte reps not worth encoding as repeat) iqueue(nqued)=ibyte nqued=nqued+1 if(ibyte.ne.iqueue(nqued-2)) then nrep=1 else c new byte is repeat nrep=nrep+1 if(nrep.eq.3) then c flush non-repeat data; change to repeat state (nrep.ge.3) c (N.B.: in this state, nqued remains at 3) nnonrep=nqued-nrep if(nnonrep.gt.0) call wb(nnonrep) do i=0,nnonrep-1 call wb(iqueue(i)) enddo irepbyte=ibyte endif endif else c currently in repeat state--does current byte continue reps? if(ibyte.ne.irepbyte) then c new byte is non-repeat--flush repeat data call wb(128+nrep) call wb(irepbyte) c place single byte in queue goto 1 endif c else, repeats continue nrep=nrep+1 endif if(nqued.lt.127 .and. nrep.lt.127) then c print *,'q: added: nr,nq:',nrep,nqued c pause ' ' return endif Entry flushq if(nrep.ge.3) then call wb(128+nrep) call wb(irepbyte) elseif (nqued.gt.0) then call wb(nqued) do i=0,nqued-1 call wb(iqueue(i)) enddo endif Entry initq nqued=0 nrep=0 c print *,'q: flush n,nq:',nrep,nqued c pause ' ' return END c***********Binary file I/O routines: write byte, string, etc.********** c (keeps track of file byte position (ipos) ) subroutine wb (i) integer*2 i,i2 integer*4 ipos common/filio/ ipos character*1 a equivalence(a,i2) i2=i write(1) a ipos=ipos+1 END subroutine wl (l) integer*4 l integer*4 ipos common/filio/ ipos write(1) l ipos=ipos+4 END subroutine ws (i) integer*2 i integer*4 ipos common/filio/ ipos write(1) i ipos=ipos+2 END subroutine wst (a) character*40 a integer*4 ipos common/filio/ ipos write (1) a(1:Len_trim(a)) ipos=ipos+Len_trim(a) END c. c. c. c*************** |