PROGRAM FSIR_LOCMAP C C THIS PROGRAM ILLUSTRATES THE READING OF A FILE IN THE BYU MERS SIR C FILE FORMAT. IT READS THE FILE, TESTS THE LOCATION TRANSFORMATIONS, C WRITES OUT A BYTE-ARRAY COPY OF THE IMAGE, AND WRITES A SIR FORMAT C OUTPUT FILE C C WRITTEN BY D.LONG: MARCH 1997 C CHARACTER*120 FNAME,FNAME2 C C IMAGE STORAGE VARIABLES C PARAMETER (MAXSIZE=4500000) ! MAXIMUM PIXELS REAL STVAL(MAXSIZE) ! READ ARRAY C LOGICAL YESNO1 C C IMAGE FILE HEADER INFORMATION C CHARACTER SENSOR*40,TITLE*80 CHARACTER TAG*40,TYPE*138,CRPROC*100,CRTIME*28 INTEGER NSX,NSY,IOPT,IDATATYPE,IOFF,ISCALE INTEGER NHEAD,NDES,NHTYPE,LDES,NIA,IPOL,IFREQHM,ISPARE1 INTEGER IREGION,ITYPE,IYEAR,ISDAY,ISMIN,IEDAY,IEMIN REAL XDEG,YDEG,ASCALE,BSCALE,A0,B0,ANODATA,VMIN,VMAX C C OPTIONAL HEADER INFO C PARAMETER (MAXI=256) CHARACTER DESCRIP*1024 DIMENSION IAOPT(MAXI) C C PROMPT USER FOR FILE NAME INPUT C WRITE(*,*) WRITE(*,15) 15 FORMAT(' Enter Input SIR File Name: ',$) READ(*,20) FNAME 20 FORMAT(A70) C 40 CONTINUE WRITE(*,*) 'SIR In= "',FNAME(1:LENGTH(FNAME)),'"' C C READ SIR FILE HEADER C MAXDES=LEN(DESCRIP) CALL READSIRHEAD(FNAME,IU,IERR,NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE, * IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX, * SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME, * MAXDES,DESCRIP,LDES,MAXI,IAOPT,NIA) C C WRITE OUT IMAGE HEADER INFORMATION C WRITE (*,*) WRITE (*,*) 'SIR File Header Information:' CALL PRINTHEAD(NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE, * IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX, * SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME, * DESCRIP,LDES,IAOPT,NIA) WRITE (*,*) C C DETERMINE IF PIXEL CENTER OR CORNER DESIRED C HALF=0.0 WRITE (*,140) 140 FORMAT(' Lat,Lon for Lower-Left corner (Y) or Center (N) of pixel? ',$) IF (.NOT.YESNO1(.TRUE.)) HALF=0.5 C C COMPUTE THE LATTITUDE IMAGE C DO IX=1,NSX DO IY=1,NSY X=IX+HALF Y=IY+HALF CALL PIX2LATLON(X,Y,THELON,THELAT,ALON,ALAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) NI=(IY-1)*NSX+IX STVAL(NI)=ALAT END DO END DO C IOFF=-91 ISCALE=360 ANODATA=-91.0 VMIN=-90.0 VMAX=90.0 C TYPE='Lat. of '//FNAME(1:LENGTH(FNAME)) FNAME2=FNAME(1:LENGTH(FNAME))//'_lat' C C WRITE SIR FORMAT FILE C NHTYPE=20 ! SET HEADER TYPE IDATATYPE=2 ! MAKE SURE OUTPUT IMAGE IS IN STANDARD I*2 FORM ITYPE=31 ! SET FILE TYPE CODE TO LATITUDE C WRITE (*,*) 'Writing "',FNAME2(1:LENGTH(FNAME2)),'"' CALL WRITESIR(FNAME2,40,IERR,NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE, * IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX, * SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME, * STVAL,DESCRIP,LDES,IAOPT,NIA) C IF (IERR.LT.0) THEN WRITE (*,*) '*** ERROR WRITING SIR FILE ***' ENDIF C C COMPUTE THE LONGITUDE IMAGE C DO IX=1,NSX DO IY=1,NSY X=IX+HALF Y=IY+HALF CALL PIX2LATLON(X,Y,THELON,THELAT,ALON,ALAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) NI=(IY-1)*NSX+IX IF (ALON.GT.180.0) ALON=ALON-360.0 IF (ALON.LT.-180.0) ALON=ALON+360.0 STVAL(NI)=ALON END DO END DO C IOFF=-181 ISCALE=180 ANODATA=-181.0 VMIN=-180.0 VMAX=180.0 C TYPE='Long. of '//FNAME(1:LENGTH(FNAME)) FNAME2=FNAME(1:LENGTH(FNAME))//'_lon' ITYPE=30 ! SET FILE TYPE CODE TO LONGITUDE C C WRITE SIR FORMAT FILE C WRITE (*,*) 'Writing "',FNAME2(1:LENGTH(FNAME2)),'"' CALL WRITESIR(FNAME2,40,IERR,NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE, * IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX, * SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME, * STVAL,DESCRIP,LDES,IAOPT,NIA) C IF (IERR.LT.0) THEN WRITE (*,*) '*** ERROR WRITING SIR FILE ***' ENDIF C 360 STOP END C C ************************************************************************** C SUBROUTINE PRINTHEAD(NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE, * IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX, * SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME, * DESCRIP,LDES,IAOPT,NIA) C C PRINTS OUT SIR FILE HEADER C CHARACTER*(*) SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,DESCRIP INTEGER IAOPT(1) C IF (NHTYPE.LT.16) WRITE (*,*) 'Old style header',nhtype C WRITE (*,*) 'Title: "',TITLE(1:LENGTH(TITLE)),'"' WRITE (*,*) 'Sensor: "',SENSOR(1:LENGTH(SENSOR)),'"' IF (NHTYPE.GT.16) THEN WRITE (*,*) 'Type: "',TYPE(1:LENGTH(TYPE)),'"' WRITE (*,*) 'Tag: "',TAG(1:LENGTH(TAG)),'"' WRITE (*,*) 'Creator:"',CRPROC(1:LENGTH(CRPROC)),'"' WRITE (*,*) 'Created:"',CRTIME(1:LENGTH(CRTIME)),'"' ENDIF IF (IOPT.EQ.-1) THEN ! IMAGE ONLY WRITE (*,*) 'Rectangular image-only form: ',ioff,iscale,iopt WRITE (*,*) ' x,y scale ',ascale,bscale WRITE (*,*) ' x,y offsets:',a0,b0 WRITE (*,*) ' x,y span: ',xdeg,ydeg ELSE IF (IOPT.EQ.0) THEN ! LAT/LONG WRITE (*,*) 'Lat/Lon Rectangular form: ',ioff,iscale,iopt WRITE (*,*) ' X and Y scales (in pixels/deg)',ascale,bscale WRITE (*,*) ' Offset values:',a0,b0 WRITE (*,*) ' Degrees: ',xdeg,ydeg ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN ! LAMBERT WRITE (*,*) 'Lambert form: ',ioff,iscale,iopt WRITE (*,*) ' X and Y scales (in km/pixel)',1./ascale,1./bscale WRITE (*,*) ' Lower Left Corner:',a0,b0 WRITE (*,*) ' Center Point: ',xdeg,ydeg IF (IOPT.EQ.2) WRITE (*,*) ' Local radius used' ELSE IF (IOPT.EQ.5) THEN ! POLAR STEREOGRAPHIC WRITE (*,*) 'Polar Sterographic form: ',ioff,iscale,iopt WRITE (*,*) ' X and Y scales (in km/pixel)',ascale,bscale WRITE (*,*) ' Lower Left Corner: (km) ',a0,b0 WRITE (*,*) ' Center Lon, Lat: ',xdeg,ydeg ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12) THEN ! EASE NORTH OR SOUTH WRITE (*,*) 'EASE polar azimuthal form: ',ioff,iscale,iopt WRITE (*,*) ' A, B scales:',ascale,bscale WRITE (*,*) ' Map origin (col/row):',a0,b0 WRITE (*,*) ' Map center (col/row):',xdeg,ydeg ELSE IF (IOPT.EQ.13) THEN ! EASE CYLINDRICAL WRITE (*,*) 'EASE cylindrical form: ',ioff,iscale,iopt WRITE (*,*) ' A, B scales:',ascale,bscale WRITE (*,*) ' Map origin (col/row):',a0,b0 WRITE (*,*) ' Map center (col/row):',xdeg,ydeg ELSE ! UNKNOWN TRANSFORMATION WRITE (*,*) '*** Unrecognized SIR file form ***',ioff,iscale,iopt WRITE (*,*) ' a,b scale ',ascale,bscale WRITE (*,*) ' a,b offsets:',a0,b0 WRITE (*,*) ' x,y span: ',xdeg,ydeg ENDIF WRITE (*,*) ' Pixels: ',nsx,' x ',nsy WRITE (*,*) ' Year: ',iyear,' Region: ',iregion WRITE (*,*) ' Start: ',isday,ismin WRITE (*,*) ' End: ',ieday,iemin WRITE (*,*) ' Type: ',itype,' Transform option: ',iopt IF (NHTYPE.GT.16) THEN WRITE (*,*) ' Frequency: ',ifreqhm*0.1,' MHz Pol: ',ipol WRITE (*,*) ' Data: ',idatatype,' Header: ',nhead,nhtype,ndes WRITE (*,*) ' Offset:',ioff,' Scale:',iscale WRITE (*,*) ' Nodata:',anodata,' Vn,Vx:',vmin,vmax IF (LDES.GT.0) WRITE (*,*) 'Description: "',DESCRIP(1:LDES),'"' IF (NIA.GT.0) THEN DO I=1,NIA WRITE (*,*) ' Aopt: ',i,' ',iaopt(i) END DO ENDIF ENDIF RETURN END C C ************************************************************************** C INTEGER FUNCTION LENGTH(LINE) C C COMPUTES THE LENGTH (NON-BLANK) OF A STRING C CHARACTER*(*) LINE N=LEN(LINE) DO I=N,1,-1 IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.CHAR(0)) THEN LENGTH=I GOTO 10 ENDIF END DO LENGTH=1 10 RETURN END C C ************************************************************************** C C STANDARD SIR IMAGE FORMAT READ/WRITE ROUTINES C C ************************************************************************** C subroutine readsirhead(fname,iu,ierr,nhead,ndes,nhtype,idatatype, * nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale, * iyear,isday,ismin,ieday,iemin,iregion,itype, * ipol,ifreqhm,ispare1,anodata,vmin,vmax, * sensor,title,type,tag,crproc,crtime, * maxdes,descrip,ldes,maxi,iaopt,nia) c c Version 2 standard SIRF image format header read routine c c written by D.G. Long 14 Mar 1997 c c inputs: c fname file name to read from c iu fortran logical unit to use c maxdes maximum number of characters in descrip c maxi maximum number of integers in iaopt c c outputs: c ierr file read error code c set to 0 for successful read of header c set to -1 for open error c set to -2 for header read error c nhead number of 512 byte blocks in header c nhtype file header type code c idatatype file data type code (0,2=i*2,1=i*1,4=i*4) def=i*2 c nsx,nsy image dimensions c xdeg..b0 projection factors c iopt projection type c ioff,iscale image storage scale factors c iyear,isday,ismin year,day,minute of start of data used in image c ieday,iemin day,minute of end of data used in image c iregion study region id number c itype image type code (=0 none or unknown) c typically: 1's digit is polarization (0=H, 1=V) c 10's digit and higher is frequency in hundreds of MHz c ispare* spare integers c sensor 20 character string: describes sensor c title 80 character string: image title c type 80 character string: description of image type c tag creator byline c crproc 100 character string: name and version of creating program c crtime 40 character string: date and time of creation c descrip variable length character string: description and annontation c ldes number of characters in descrip c iaopt optional integers c by convention, first value iaopt is a code telling how to interpret c the rest of the array if nia>0 c nia number of integers c c SIR files have headers which are multiples of 512 bytes followed by c image data stored as 2 byte integers zero padded to be a multiple c of 512 bytes. c character*(*) fname integer iu,nhead,nhtype,idatatype integer nsx,nsy,iopt real xdeg,ydeg,ascale,bscale,a0,b0 integer ioff,iscale,iyear,isday,ismin,ieday,iemin integer iregion,itype,ipol,ifreqhm,ispare1 real anodata,vmin,vmax character*(*) sensor,title,type,tag,crproc,crtime,descrip integer iaopt(*) c c working variables c integer*2 i2(2) real f2 equivalence (i2,f2) ! used in i*2 to real mapping integer*2 temp(256) character sens1*41,tit1*81,type1*139 character tag1*41,crp1*101,crt1*29 c c note: the bswap and vms flags should be set in the block data c logical bswap,vms common /vaxorunix/bswap,vms save /vaxorunix/ c c header mapping (see code for scaling of <= variables) c note that character strings do not have to be null terminated c c temp(1) = nsx ! pixels in x direction c temp(2) = nsy ! pixels in y direction c temp(3) <= xdeg ! span of deg in x c temp(4) <= ydeg ! span of deg in y c temp(5) = nhtype ! header type (old<15) c temp(6) <= ascale ! x scaling c temp(7) <= bscale ! y scaling c temp(8) <= a0 ! x (lon) origin c note: longitudes should be in the range -180 to 180 c temp(9) <= b0 ! y (lat) origin c temp(10) = ioff ! offset to be added to scaled val c temp(11) = iscale ! scale factor ival=(val-ioff)/iscale c temp(12) = iyear ! year for data used c temp(13) = isday ! starting JD c temp(14) = ismin ! time of day for first data (in min) c temp(15) = ieday ! ending JD c temp(16) = iemin ! time of day for last data (in min) c temp(17) = iopt ! projection type c ! -1 = no projection, image only c ! 0 = rectalinear lat/lon c ! 1 = lambert equal area c ! 2 = lambert equal area (local rad) c ! 5 = polar stereographic c ! 11 = EASE north equal area grid c ! 12 = EASE south equal area grid c ! 13 = EASE cylindrical grid c temp(18) = iregion ! region id code c temp(19) = itype ! image type code c ! standard values: 0=unknown or n/a c ! 1 = scatterometer A (dB) c ! 2 = scatterometer B (dB/deg) c ! 3 = radiometer Tb (K) c ! 9 = topography (m) c temp(20)-temp(39) 40 chars of sensor c temp(40) = 0 c temp(41) = nhead ! number of 512 byte header blocks c temp(42) = ndes ! number of 512 byte blocks description c temp(43) = ldes ! number of bytes of description c temp(44) = nia ! number of optional integers c temp(45) = ipol ! polarization (0=n/a,1=H,2=V) c temp(46) = ifreqhm ! frequency in 100's MHz (0 if n/a) c temp(47) = ispare1 ! spare c temp(48) = idatatype ! data type code 0,2=i*2,1=i*1,4=f c c the value of idata type determines how data is stored and how c anodata, vmin, and vmax are stored. c c if idatatype = 1 data is stored as bytes and minv=128 c if idatatype = 2 data is stored as 2 byte integers and minv=32766 c if idatatype = 4 data is stored as IEEE floating point c c if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers c in temp(49)..temp(51) minv, ioff and iscal used to convert c integers or bytes into floating point values c nodata, vmin, and vmax must be representable with ioff and iscale c temp(*) = (value-ioff)*iscale-minv c value = float(temp(*)+minv)/float(iscale)+ioff c idatatype=2 is considered the SIR standard format c c if idatatype = f anodata,vmin,vmax are stored as floating points c in temp(52)..temp(57) and minv, ioff and iscale are ignored here c and when reading the file. c floating point numbers are NOT standard across platforms and c is not recommended c c temp(49) <= anodata ! value representing no data c temp(50) <= vmin ! minimum useful value from creator prg c temp(51) <= vmax ! maximum useful value from creator prg c temp(52,53) = anodata ! IEEE floating value of no data c temp(54,55) = vmin ! IEEE floating minimum useful value c temp(56,57) = vmax ! IEEE floating maximum useful value c c temp(58)-temp(126) 138 chars of type c temp(127) = 0 c temp(128) = 0 c temp(129)-temp(168) 80 chars of title c temp(169) = 0 c temp(170)-temp(189) 40 chars of tag c temp(190) = 0 c temp(191)-temp(240) 100 chars of crproc c temp(241) = 0 c temp(242)-temp(255) 28 chars of crtime c temp(256) = 0 c c optional header blocks: c c ndes header blocks of 512 bytes: chars of description c nhead-ndes-1 header blocks of 512 bytes: values of iaopt c by convention, first value iaopt is a code telling how to interpret c the rest of the array if nia>0 c c remainder of file is image data in a multiple of 512 byte blocks c two byte integer and byte scaling is c intval = (fvalue-ioff)*iscale-minv c fvalue = float(intval+minv)/float(iscale)+ioff c c open input file c if (vms) then open(unit=iu,name=fname,status='old',form='formatted', * recl=512,err=99) !, c * readonly,recordtype='fixed',carriagecontrol='none') else open(unit=iu,name=fname,status='old',form='unformatted', * recl=512,access='direct',err=99) endif c c read first header block c if (vms) then read (iu,414) temp 414 format(256a2) call swapbuf(temp,256) else read (iu,rec=1) temp if (bswap) call swapbuf(temp,256) endif c c decode first header block c nsx = temp(1) ! dimension nsy = temp(2) iopt = temp(17) ! transformation option nhtype = temp(5) ! header type if (nhtype.lt.20) nhtype=1 c c get default transformation decoding c xdeg = 0.01 * temp(3) ! center point or span ydeg = 0.01 * temp(4) ascale = temp(6) * 0.001 ! scale factor bscale = temp(7) * 0.001 a0 = temp(8) * 0.01 ! origin b0 = temp(9) * 0.01 c c get special cases which depend on transformation option c if (iopt.eq.-1) then ! image only xdeg = 0.1 * temp(3) ydeg = 0.1 * temp(4) else if (iopt.eq.0) then ! rectalinear xdeg = 0.01 * temp(3)+100.0 else if (iopt.eq.1.or.iopt.eq.2) then ! lambert ascale=1000.0/float(temp(6)) bscale=1000.0/float(temp(7)) a0=temp(8) b0=temp(9) else if (iopt.eq.5) then ! polar stereographic ascale = temp(6) * 0.01 bscale = temp(7) * 0.01 a0 = temp(8) b0 = temp(9) xdeg = 0.01 * temp(3)+100.0 else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE grid ascale = 2.0d0*(0.001*float(temp(6)))*6371.228d0/ * 25.067525d0 bscale = 2.0d0*(0.001*float(temp(7)))*25.067525d0 xdeg = 0.1 * temp(3) ydeg = 0.1 * temp(4) a0 = 0.1 * temp(8) b0 = 0.1 * temp(9) else write (*,*)'*** Unrecognized SIR option in READSIRHEAD ***' endif c c decode scaling, data time, region, and type information c ioff = temp(10) iscale = temp(11) if (iscale.eq.0) iscale=1 iyear = temp(12) isday = temp(13) ismin = temp(14) ieday = temp(15) iemin = temp(16) iregion = temp(18) itype = temp(19) nhead = temp(41) ! total number of header blocks if (nhead.eq.0) nhead=1 if (nhtype.eq.1) nhead=1 ndes = temp(42) ! number of description blocks ldes = temp(43) nia = temp(44) ipol = temp(45) ifreqhm = temp(46) ispare1 = temp(47) idatatype = temp(48) ! data type code if (idatatype.eq.0) idatatype=2 soff=32767.0/float(iscale) if (idatatype.eq.1) soff=128.0/float(iscale) anodata = temp(49)/float(iscale)+ioff+soff vmin = temp(50)/float(iscale)+ioff+soff vmax = temp(51)/float(iscale)+ioff+soff c if (idatatype.eq.4) then i2(1)=temp(52) i2(2)=temp(53) anodata=f2 i2(1)=temp(54) i2(2)=temp(55) vmin=f2 i2(1)=temp(56) i2(2)=temp(57) vmax=f2 endif c c get character variables c do i=1,20 j=(i-1)*2+1 sens1(j:j)=char(mod(temp(i+19),256)) sens1(j+1:j+1)=char(temp(i+19)/256) end do sens1(41:41)=char(0) sensor=sens1 do i=1,40 j=(i-1)*2+1 tit1(j:j)=char(mod(temp(i+128),256)) tit1(j+1:j+1)=char(temp(i+128)/256) end do tit1(81:81)=char(0) title=tit1 if (nhtype.gt.1) then do i=1,69 j=(i-1)*2+1 type1(j:j)=char(mod(temp(i+57),256)) type1(j+1:j+1)=char(temp(i+57)/256) end do type1(139:139)=char(0) do i=1,20 j=(i-1)*2+1 tag1(j:j)=char(mod(temp(i+169),256)) tag1(j+1:j+1)=char(temp(i+169)/256) end do tag1(41:41)=char(0) do i=1,50 j=(i-1)*2+1 crp1(j:j)=char(mod(temp(i+190),256)) crp1(j+1:j+1)=char(temp(i+190)/256) end do crp1(101:101)=char(0) do i=1,14 j=(i-1)*2+1 crt1(j:j)=char(mod(temp(i+241),256)) crt1(j+1:j+1)=char(temp(i+241)/256) end do crt1(29:29)=char(0) else type1=' ' tag1=' ' crp1=' ' crt1=' ' endif type=type1 tag=tag1 crproc=crp1 crtime=crt1 descrip=' ' c c get extra storage information c if (nhead.gt.1) then if (ndes.gt.0) then c c read next header blocks c icnt=0 do irec=2,ndes+1 if (vms) then read (iu,414) temp call swapbuf(temp,256) else read (iu,rec=irec) temp if (bswap) call swapbuf(temp,256) endif do j=1,512 icnt=icnt+1 if (icnt.le.ldes.and.icnt.le.maxdes) then i=(j-1)/2+1 iv=temp(i) if (mod(j,2).eq.0) then iv=iv/256 else iv=mod(iv,256) endif descrip(icnt:icnt)=char(iv) endif end do end do if (maxdes.lt.ldes) ldes=maxdes endif if (nhead-ndes-1.gt.0) then icnt=0 do irec=ndes+2,nhead if (vms) then read (iu,414) temp call swapbuf(temp,256) else read (iu,rec=irec) temp if (bswap) call swapbuf(temp,256) endif do i=1,256 icnt=icnt+1 if (icnt.le.nia.and.icnt.le.maxi) * iaopt(icnt) = temp(i) end do end do if (nia.gt.maxi) nia=maxi endif endif return c 300 continue write (*,398) 398 format(' *** ERROR reading header of input SIR file ***') c c close file on header read error c close(unit=iu) c ierr=-2 return 99 continue write (*,198) 198 format(' *** ERROR opening input SIR file ***') ierr=-1 return end c c subroutine readsirf(iu,ierr,nhead,nhtype,idatatype,stval,nsx,nsy, * ioff,iscale,smin,smax,ncnt,anodata,vmin,vmax) c c read image data c c inputs: c iu fortran logical unit of previously opened file c ierr error code from header read c nhead number of 512 byte blocks in header c nhtype header type c nsx,nsy image dimensions c ioff,iscale image storage scale factors c anodata specify nodata value c vmin,vmax specify min, max valid range c c outputs: c ierr file read error code c set to 0 for successful read of file c set to -3 for data read error c stval output image array c smin,max min and max of image within valid range c valid range = [vmin..vmax] and .ne. anodata c ncnt count of pixels in valid range c integer iu,nhead,nhtype real stval(*) integer nsx,nsy,ncnt integer ioff,iscale real anodata,vmin,vmax,smin,smax c c working variables c byte tempb(512) integer*2 temp(256) real tempf(128) integer ntot,jj,irec real s,soff c c note: the bswap and vms flags should be set in the block data c logical bswap,vms common /vaxorunix/bswap,vms save /vaxorunix/ c if (ierr.lt.0) return ! return if file open previously failed c ntot = nsx*nsy if (iscale.eq.0) iscale=1 s=1./float(iscale) soff=32767.0/float(iscale) nrec=256 if (idatatype.eq.1) then soff=128.0/float(iscale) nrec=512 endif if (idatatype.eq.4) nrec=128 c jj=0 irec=nhead smin=1.e25 smax=-1.e25 ncnt=0 200 continue irec=irec+1 if (vms) then if (idatatype.eq.1) then read (iu,413,end=299) tempb 413 format(512a1) else if (idatatype.eq.4) then read (iu,412,end=299) tempf 412 format(128a4) else read (iu,414,end=299) temp 414 format(256a2) call swapbuf(temp,256) endif else if (idatatype.eq.1) then read(iu,rec=irec,err=298) tempb else if (idatatype.eq.4) then read(iu,rec=irec,err=298) tempf else read(iu,rec=irec,err=298) temp if (bswap) call swapbuf(temp,256) endif endif do i=1,nrec jj=jj+1 if (idatatype.eq.1) then k=tempb(i) stval(jj)=s*k+soff+ioff else if (idatatype.eq.4) then stval(jj)=tempf(i) else stval(jj)=s*temp(i)+soff+ioff endif if (stval(jj).ne.anodata.and. * stval(jj).ge.vmin.and. * stval(jj).le.vmax) then smin=min(smin,stval(jj)) smax=max(smax,stval(jj)) ncnt=ncnt+1 endif if (jj.ge.ntot) goto 299 end do goto 200 298 continue write (*,98) 98 format(' *** ERROR reading input SIR file ***') c c close file on read error c close(unit=iu) ierr=-3 return 299 continue c c close file on successful completion c close(unit=iu) ierr=0 return end c c ************************************************************************* c c standard sir write routine c subroutine writesir(fname,iu,ierr,nhead,ndes,nhtype,idatatype, * nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale, * iyear,isday,ismin,ieday,iemin,iregion,itype, * ipol,ifreqhm,ispare1,anodata,vmin,vmax, * sensor,title,type,tag,crproc,crtime, * stval,descrip,ldes,iaopt,nia) c c Version 2 standard SIRF image format write routine c c written by D.G. Long 15 Mar 1997 c c inputs: c fname file name to read from c iu fortran logical unit to use c idatatype file data type code (0,2=i*2,1=i*1,4=i*4) def=i*2 c nsx,nsy image dimensions c xdeg..b0 projection factors c iopt projection type c ioff,iscale image storage scale factors c iyear,isday,ismin year,day,minute of start of data used in image c ieday,iemin day,minute of end of data used in image c iregion study region id number c itype image type code (=0 none or unknown) c ispare* spare integers c sensor 20 character string: describes sensor c title 80 character string: image title c type 80 character string: description of image type c tag creator byline c crproc 100 character string: name and version of creating program c crtime 40 character string: date and time of creation c stval image array c descrip variable length character string: description and annontation c ldes number of characters in descrip c iaopt optional integers c by convention, first value iaopt is a code telling how to interpret c the rest of the array if nia>0 c nia number of integers c c outputs: c ierr file write error code c set to 0 for successful write c set to -1 for open error c set to -2 for write error c nhead number of 512 byte blocks in header c nhtype file header type code c c SIR files have headers which are multiples of 512 bytes followed by c image data stored as 2 byte integers zero padded to be a multiple c of 512 bytes. c character*(*) fname integer iu,nhead,nhtype,idatatype integer nsx,nsy,iopt real xdeg,ydeg,ascale,bscale,a0,b0 integer ioff,iscale,iyear,isday,ismin,ieday,iemin integer iregion,itype,ipol,ifreqhm,ispare1 real anodata,vmin,vmax real stval(*) character*(*) sensor,title,type,tag,crproc,crtime,descrip integer iaopt(*) c c working variables c byte tempb(512) integer*2 temp(256) real tempf(128) equivalence (temp,tempb) ! used only to save space equivalence (temp,tempf) ! used only to save space integer*2 i2(2) real f2 equivalence (f2,i2) ! used in i*2 to float mapping character sens1*41,tit1*81,type1*139 character tag1*41,crp1*101,crt1*29 c c note: the bswap and vms flags should be set in the block data c logical bswap,vms common /vaxorunix/bswap,vms save /vaxorunix/ c c header mapping (see code for scaling of <= variables) c note that character strings do not have to be null terminated c c temp(1) = nsx ! pixels in x direction c temp(2) = nsy ! pixels in y direction c temp(3) <= xdeg ! span of deg in x c temp(4) <= ydeg ! span of deg in y c temp(5) = nhtype ! header type (old<15) c temp(6) <= ascale ! x scaling c temp(7) <= bscale ! y scaling c temp(8) <= a0 ! x (lon) origin c note: longitudes should be in the range -180 to 180 c temp(9) <= b0 ! y (lat) origin c temp(10) = ioff ! offset to be added to scaled val c temp(11) = iscale ! scale factor ival=(val-ioff)/iscale c temp(12) = iyear ! year for data used c temp(13) = isday ! starting JD c temp(14) = ismin ! time of day for first data (in min) c temp(15) = ieday ! ending JD c temp(16) = iemin ! time of day for last data (in min) c temp(17) = iopt ! projection type c ! -1 = no projection, image only c ! 0 = rectalinear lat/lon c ! 1 = lambert equal area c ! 2 = lambert equal area (local rad) c ! 5 = polar stereographic c ! 11 = EASE north equal area grid c ! 12 = EASE south equal area grid c ! 13 = EASE cylindrical grid c temp(18) = iregion ! region id code c temp(19) = itype ! image type code c ! standard values: 0=unknown or n/a c ! 1 = scatterometer A (dB) c ! 2 = scatterometer B (dB/deg) c ! 3 = radiometer Tb (K) c ! 9 = topography (m) c temp(20)-temp(39) 40 chars of sensor c temp(40) = 0 c temp(41) = nhead ! number of 512 byte header blocks c temp(42) = ndes ! number of 512 byte blocks description c temp(43) = ldes ! number of bytes of description c temp(44) = nia ! number of optional integers c temp(45) = ipol ! polarization (0=n/a,1=H,2=V) c temp(46) = ifreqhm ! frequency in 100's MHz (0 if n/a) c temp(47) = ispare1 ! spare c temp(48) = idatatype ! data type code 0,2=i*2,1=i*1,4=i*4 c c the value of idata type determines how data is stored and how c anodata, vmin, and vmax are stored. c c if idatatype = 1 data is stored as bytes and minv=128 c if idatatype = 2 data is stored as 2 byte integers and minv=32766 c if idatatype = 4 data is stored as IEEE floating point c c if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers c in temp(49)..temp(51) minv, ioff and iscal used to convert c integers or bytes into floating point values c nodata, vmin, and vmax must be representable with ioff and iscale c temp(*) = (value-ioff)*iscale-minv c value = float(temp(*)+minv)/float(iscale)+ioff c idatatype=2 is considered the SIR standard format c c if idatatype = f anodata,vmin,vmax are stored as floating points c in temp(42)..temp(57) and minv, ioff and iscale are ignored here c and when reading the file. c floating point numbers are NOT standard across platforms and c is not recommended c c temp(49) <= anodata ! value representing no data c temp(50) <= vmin ! minimum useful value from creator prg c temp(51) <= vmax ! maximum useful value from creator prg c temp(52,53) = anodata ! IEEE floating value of no data c temp(54,55) = vmin ! IEEE floating minimum useful value c temp(56,57) = vmax ! IEEE floating maximum useful value c c temp(58)-temp(126) 150 chars of type c temp(127) = 0 c temp(128) = 0 c temp(129)-temp(168) 80 chars of title c temp(169) = 0 c temp(170)-temp(189) 40 chars of tag c temp(190) = 0 c temp(191)-temp(240) 100 chars of crproc c temp(241) = 0 c temp(242)-temp(255) 28 chars of crtime c temp(256) = 0 c c optional header blocks: c c ndes header blocks of 512 bytes: chars of description c nhead-ndes-1 header blocks of 512 bytes: values of iaopt c by convention, first value iaopt is a code telling how to interpret c the rest of the array if nia>0 c c remainder of file is image data in a multiple of 512 byte blocks c two byte integer scaling is c intval = (fvalue-ioff)*iscale-32767 c fvalue = float(intval+32767)/float(iscale)+ioff c c make first header block c ierr = 0 temp(1) = nsx ! image dimensions temp(2) = nsy temp(5) = 20 ! version 2.0 header c c scale image transformation variables based on transform option c if (iopt.eq.-1) then ! image only temp(3) = nint(xdeg * 10.0) temp(4) = nint(ydeg * 10.0) temp(6) = nint(1000. * ascale) temp(7) = nint(1000. * bscale) temp(8) = nint(100. * a0) temp(9) = nint(100. * b0) else if (iopt.eq.0) then ! rectalinear lat/lon temp(3) = nint((xdeg-100.0) * 100.0) ! lon +/- 180.0 deg temp(4) = nint(ydeg * 100.0) temp(6) = nint(1000.*ascale) temp(7) = nint(1000.*bscale) temp(8) = nint(100.*a0) temp(9) = nint(100.*b0) else if (iopt.eq.1.or.iopt.eq.2) then ! lambert temp(3) = nint(xdeg * 100.0) ! lon +/- 180.0 deg temp(4) = nint(ydeg * 100.0) temp(6) = nint(1000./ascale) temp(7) = nint(1000./bscale) temp(8) = nint(a0) temp(9) = nint(b0) else if (iopt.eq.5) then ! polar stereographic temp(3) = nint((xdeg-100.0) * 100.0) ! lon +/- 180.0 deg temp(4) = nint(ydeg * 100.0) temp(6) = nint(100.*ascale) temp(7) = nint(100.*bscale) temp(8) = nint(a0) temp(9) = nint(b0) else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE grid temp(3) = nint(xdeg * 10.0) temp(4) = nint(ydeg * 10.0) temp(6) = nint(1000.*anint(10.*ascale*25.067525/ * 6371.228)*0.05) temp(7) = nint(1000.*anint(10.*bscale/25.067525)*0.05) temp(8) = nint(10. * a0) temp(9) = nint(10. * b0) else temp(3) = nint(xdeg * 100.0) temp(4) = nint(ydeg * 100.0) temp(6) = nint(1000. * ascale) temp(7) = nint(1000. * bscale) temp(8) = nint(100. * a0) temp(9) = nint(100. * b0) write (*,*) '*** Unrecognized SIR option in WRITESIR ***' endif c c store scale factors, time, option, region, and image type c temp(10) = ioff temp(11) = iscale if (iscale.eq.0) temp(11)=1 temp(12) = iyear temp(13) = isday temp(14) = ismin temp(15) = ieday temp(16) = iemin temp(17) = iopt temp(18) = iregion temp(19) = itype c nhead=1 ndes=0 if (ldes.gt.0) then ndes=ldes/512 if (mod(ldes,512).gt.0) ndes=ndes+1 write (*,*) 'extra descriptor',ldes,ndes nhead=nhead+ndes endif if (nia.gt.0) then write (*,*) 'extra integers',nia nhead=nhead+nia/256 if (mod(nia,256).gt.0) nhead=nhead+1 endif c temp(41) = nhead temp(42) = ndes temp(43) = ldes temp(44) = nia temp(45) = ipol temp(46) = ifreqhm temp(47) = ispare1 temp(48) = idatatype c if (idatatype.eq.1) then maxv=127 minv=128 else maxv=32767 minv=32767 endif c if (idatatype.ne.4) then k = nint(iscale*(anodata-ioff))-minv if (k.gt.maxv) then k=maxv write (*,*) '*** WRITESIR: overflow in "anodata"',anodata,ioff,iscale endif if (k.lt.-minv) then k=-minv type *,'anodata',anodata,vmin,ioff,iscale,k,minv write (*,*) '*** WRITESIR: underflow in "anodata"',anodata,ioff,iscale endif temp(49) = k k = nint(iscale*(vmin-ioff))-minv if (k.gt.maxv) then k=maxv write (*,*) '*** WRITESIR: overflow in "vmin"',vmin,ioff,iscale endif if (k.lt.-minv) then k=-minv write (*,*) '*** WRITESIR: underflow in "vmin"',vmin,ioff,iscale endif temp(50) = k k = nint(iscale*(vmax-ioff))-minv if (k.gt.maxv) then k=maxv write (*,*) '*** WRITESIR: overflow in "vmax"',vmax,ioff,iscale endif if (k.lt.-minv) then k=-minv write (*,*) '*** WRITESIR: underflow in "vmax"',vmax,ioff,iscale endif temp(51) = k endif c f2=anodata temp(52)=i2(1) temp(53)=i2(2) f2=vmin temp(54)=i2(1) temp(55)=i2(2) f2=vmax temp(56)=i2(1) temp(57)=i2(2) c c write character variables c sens1=' ' tit1=' ' type1=' ' tag1=' ' crp1=' ' crt1=' ' sens1=sensor do i=1,20 j=(i-1)*2+1 temp(i+19)=ichar(sens1(j:j))+256*ichar(sens1(j+1:j+1)) end do temp(40) = 0 type1=type do i=1,69 j=(i-1)*2+1 temp(i+57)=ichar(type1(j:j))+256*ichar(type1(j+1:j+1)) end do temp(127) = 0 temp(128) = 0 tit1=title do i=1,40 j=(i-1)*2+1 temp(i+128)=ichar(tit1(j:j))+256*ichar(tit1(j+1:j+1)) end do temp(169) = 0 tag1=tag do i=1,20 j=(i-1)*2+1 temp(i+169)=ichar(tag1(j:j))+256*ichar(tag1(j+1:j+1)) end do temp(190) = 0 crp1=crproc do i=1,50 j=(i-1)*2+1 temp(i+190)=ichar(crp1(j:j))+256*ichar(crp1(j+1:j+1)) end do temp(241) = 0 crt1=crtime do i=1,14 j=(i-1)*2+1 temp(i+241)=ichar(crt1(j:j))+256*ichar(crt1(j+1:j+1)) end do temp(256) = 0 c c open output file c if (vms) then open(unit=iu,name=fname,status='new',form='formatted',err=999, * recl=512) c * ,recordtype='fixed',carriagecontrol='none') else open(unit=iu,name=fname,status='unknown',form='unformatted',err=999, * recl=512,access='direct') endif c c write first header block c irec=1 if (bswap) call swapbuf(temp,256) if (vms) then write(iu,414,err=1999) temp else write(iu,rec=1,err=1999) temp endif c c write additional header records if needed c if (ndes.gt.0) then icnt=0 do ii=1,ndes do i=1,256 k=0 icnt=icnt+1 if (icnt.le.ldes) then k=ichar(descrip(icnt:icnt)) icnt=icnt+1 if (icnt.le.ldes) k=k+256*ichar(descrip(icnt:icnt)) endif temp(i)=k end do if (bswap) call swapbuf(temp,256) irec=irec+1 if (vms) then write(iu,414,err=1999) temp else write(iu,rec=irec,err=1999) temp endif end do endif c c extra parameters c if (nhead-ndes-1.gt.0) then icnt=0 do ii=nhead-ndes+1,nhead do i=1,256 icnt=icnt+1 if (icnt.le.nia) then temp(i)=iaopt(icnt) else temp(i)=0 endif end do if (bswap) call swapbuf(temp,256) irec=irec+1 if (vms) then write(iu,414,err=1999) temp else write(iu,rec=irec,err=1999) temp endif end do endif c c write image data c iucnt=0 iocnt=0 ntot = nsx*nsy nrec=256 if (idatatype.eq.1) nrec=512 if (idatatype.eq.4) nrec=128 do j = 0,ntot,nrec imax=min(nrec,ntot-j) if (imax.gt.0) then do i=1,imax if (idatatype.eq.4) then tempf(i)=stval(j+i) else k = nint(iscale*(stval(j+i)-ioff))-minv if (k.gt.maxv) then k=maxv iocnt=iocnt+1 endif if (k.lt.-minv) then k=-minv iucnt=iucnt+1 endif if (idatatype.eq.1) then tempb(i) = k else temp(i) = k endif endif end do if (imax.ne.nrec) then do i=imax+1,nrec if (idatatype.eq.1) then tempb(i) = -minv else if (idatatype.eq.4) then tempf(i) = anodata else temp(i) = -minv endif end do endif irec=irec+1 if (idatatype.eq.1) then ! byte size if (vms) then write(iu,413,err=1999) tempb 413 format(512a1) else write(iu,rec=irec,err=1999) tempb endif else if (idatatype.eq.4) then ! 4 byte float size if (vms) then write(iu,412,err=1999) tempf 412 format(128a4) else write(iu,rec=irec,err=1999) tempf endif else ! integer*2 size if (bswap) call swapbuf(temp,256) if (vms) then write(iu,414,err=1999) temp 414 format(256a2) else write(iu,rec=irec,err=1999) temp endif endif endif end do c c close file c close(unit=iu) c if (iucnt.gt.0) write (*,50) iucnt,ntot 50 format(' *** WRITESIR Underflow count: ',I,' of ',I,' ***') if (iocnt.gt.0) write (*,60) iocnt,ntot 60 format(' *** WRITESIR Overflow count: ',I,' of ',I,' ***') return 999 continue write (*,5000) 5000 format(' *** ERROR OPENING OUTPUT FILE ***') ierr=-1 return 1999 continue write (*,5050) 5050 format(' *** ERROR WRITING OUTPUT FILE ***') ierr=-2 return end c c subroutine swapbuf(t,n) c c swap byte order of 2 byte integers c logical*1 t(1),s do i=1,n j=(i-1)*2+1 s=t(j) t(j)=t(j+1) t(j+1)=s end do return end c c ************************************************************************* c c this block data defines the file access and byte order logicals c in their common block. set the parameters according to the machine c code is execute. c BLOCK DATA LOGICAL BSWAP,VMS COMMON /VAXORUNIX/BSWAP,VMS C DATA VMS/.TRUE./ ! VAXVMS FILE ACCESS & BYTE ORDER DATA VMS/.FALSE./ ! UNIX FILE ACCESS DATA BSWAP/.FALSE./ ! HP, SUN (big-endian) C DATA BSWAP/.TRUE./ ! VAX, PC (little-endian) SAVE /VAXORUNIX/ END C C *************************************************************************** C C IMAGE PIXEL TO LAT,LONG AND LAT,LONG TO IMAGE PIXEL C COORDINATE TRANSFORMATION ROUTINES. NOTE THAT SIR FILE IMAGE C HEADER INFORMATION IS PASSED VIA THE COMMON BLOCK /IMAGEINFO/ C INTO THE HIGH LEVEL ROUTINES. C C *************************************************************************** C SUBROUTINE PIX2LATLON(X,Y,THELON,THELAT,ALON,ALAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) C C GIVEN AN IMAGE PIXEL LOCATION (X,Y) (AS FLOATING POINT NUMBERS) C COMPUTES THE X,Y TRANSFORMATION COORDINATES (THELON,THELAT) AND C THE LAT,LON COORDINATES (ALON,ALAT). THE LAT,LON RETURNED C CORRESPONDS TO THE LOWER-LEFT CORNER OF THE PIXEL. IF LAT,LON C OF PIXEL CENTER IS DESIRED USE (X+0.5,Y+0.5) WHERE X,Y ARE INTEGER C VALUED. C C NOTE: GIVEN INTEGER PIXEL INDICES, SET X=FLOAT(IX) Y=FLOAT(IY). C -- SEE LATLON2PIX. WHILE ROUTINE WILL ATTEMPT TO CONVERT ANY (X,Y) C VALUES, ONLY (X,Y) VALUES WITH 1 <= X <= NSX+1 AND 1 <= Y <= NSY+1 C ARE CONTAINED WITHIN IMAGE PIXELS. C C WRITTEN BY: DGL 4 FEB 1995 C REAL X,Y,THELON,THELAT,ALON,ALAT,XDEG,YDEG,ASCALE,BSCALE,A0,B0 INTEGER IOPT C IF (IOPT.EQ.-1) THEN ! IMAGE ONLY (CAN'T TRANSFORM!) THELON=(X-1.0)/ASCALE+A0 THELAT=(Y-1.0)/BSCALE+B0 ALON=THELON ALAT=THELAT ELSE IF (IOPT.EQ.0) THEN ! LAT/LONG (ASCALE AND BSCALE > 0) THELON=(X-1.0)/ASCALE+A0 THELAT=(Y-1.0)/BSCALE+B0 ALON=THELON ALAT=THELAT ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN ! LAMBERT (ASCALE AND BSCALE < 0) THELON=(X-1.0)/ASCALE+A0 THELAT=(Y-1.0)/BSCALE+B0 CALL ILAMBERT1(ALAT,ALON,THELON,THELAT,YDEG,XDEG,IOPT) ELSE IF (IOPT.EQ.5) THEN ! POLAR STEREOGRAPHIC THELON=(X-1.0)*ASCALE+A0 THELAT=(Y-1.0)*BSCALE+B0 CALL IPOLSTER(ALON,ALAT,THELON,THELAT,XDEG,YDEG) ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12.OR.IOPT.EQ.13) THEN ! EASE GRID THELON=X-XDEG-1.0-(XDEG+A0) THELAT=Y-YDEG-1.0-(YDEG+B0) CALL IEASEGRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE) ELSE ! UNKNOWN TRANSFORM ENDIF RETURN END C C SUBROUTINE LATLON2PIX(ALON,ALAT,X,Y,THELON,THELAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) C C CONVERT A LAT,LON COORDINATE (ALON,ALAT) TO AN IMAGE PIXEL LOCATION C (X,Y) (IN FLOATING POINT) AND IMAGE X,Y COORDINATE (THELON,THELAT). C TO COMPUTE INTEGER PIXEL INDICES (IX,IY): CHECK TO INSURE C 1 <= X < NSX+1 AND 1 <= X < NSX+1 THEN IX=IFIX(X) IY=IFIX(Y) C C WRITTEN BY: DGL 4 FEB 1995 C REAL X,Y,THELON,THELAT,ALON,ALAT,XDEG,YDEG,ASCALE,BSCALE,A0,B0 INTEGER IOPT C IF (IOPT.EQ.-1) THEN ! IMAGE ONLY (CAN'T TRANSFORM!) X=ASCALE*(ALON-A0)+1.0 Y=BSCALE*(ALAT-B0)+1.0 THELON=ALON THELAT=ALAT ELSE IF (IOPT.EQ.0) THEN ! LAT/LONG (ASCALE AND BSCALE > 0) THELON=ALON THELAT=ALAT X=ASCALE*(THELON-A0)+1.0 Y=BSCALE*(THELAT-B0)+1.0 ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN ! LAMBERT (ASCALE AND BSCALE < 0) CALL LAMBERT1(ALAT,ALON,THELON,THELAT,YDEG,XDEG,IOPT) X=(THELON-A0)*ASCALE+1.0 Y=(THELAT-B0)*BSCALE+1.0 ELSE IF (IOPT.EQ.5) THEN ! POLAR STEREOGRAPHIC CALL POLSTER(ALON,ALAT,THELON,THELAT,XDEG,YDEG) X=(THELON-A0)/ASCALE+1.0 Y=(THELAT-B0)/BSCALE+1.0 ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12.OR.IOPT.EQ.13) THEN ! EASE GRID CALL EASEGRID(IOPT,ALAT,ALON,THELON,THELAT,ASCALE) THELON=THELON+XDEG THELAT=THELAT+YDEG X=THELON+1.0-(XDEG+A0) Y=THELAT+1.0-(YDEG+B0) ELSE ! UNKNOWN TRANSFORMATION TYPE ENDIF C RETURN END C C SUBROUTINE F2IPIX(X,Y,IX,IY,NSX,NSY) C C QUANTIZES F FLOATING POINT PIXEL VALUE TO THE "ACTUAL" INTEGER C PIXEL VALUE (SEE LATLON2PIX) C IF (X.GE.1.0.AND.X.LT.NSX+1) THEN IX=IFIX(X+0.0001) ! ALLOW A SMALL AMOUNT OF ROUNDING ELSE IX=0 ENDIF IF (Y.GE.1.0.AND.Y.LT.NSY+1) THEN IY=IFIX(Y+0.0001) ! ALLOW A SMALL AMOUNT OF ROUNDING ELSE IY=0 ENDIF RETURN END C C ********************************************************************* C C PROJECTION TRANSFORMATION ROUTINES C C ********************************************************************* C SUBROUTINE LAMBERT1(LAT,LON,X,Y,ORGLAT,ORGLON,IOPT) C C COMPUTES THE TRANSFORMATION FROM LAT/LON TO X/Y FOR THE C LAMBERT AZIMUTHAL EQUAL-AREA PROJECTION C C SEE "MAP PROJECTIONS USED BY THE U.S. GEOLOGICAL SURVEY" C GEOLOGICAL SURVEY BULLETIN 1532, PGS 157-173 C C FOR THIS ROUTINE, A SPHERICAL EARTH IS ASSUMED FOR THE PROJECTION. C THE ERROR WILL BE SMALL FOR SMALL-SCALE MAPS. C FOR IOPT=1 A FIXED, NOMINAL EARTH RADIUS IS USED. C FOR IOPT=2 THE LOCAL RADIUS OF THE EARTH IS USED BASED ON C THE 1972 WGS ELLIPSOID MODEL (BULLETIN PG 15). C C WRITTEN BY: DGL 7 FEB 1992 C MODIFIED: DGL 4 MAR 1995 C + INCLUDED RADIUS OF THE EARTH COMPUTED AT LOCAL REFERENCE POINT C C INPUTS: C LAT (R): LATITUDE +90 TO -90 DEG WITH NORTH POSITIVE C LON (R): LONGITUDE 0 TO +360 DEG WITH EAST POSITIVE C OR -180 TO +180 WITH EAST MORE POSITIVE C ORGLAT (R): ORIGIN PARALLEL +90 TO -90 DEG WITH NORTH POSITIVE C ORGLON (R): CENTRAL MERIDIAN (LONGITUDE) 0 TO +360 DEG C OR -180 TO +180 WITH EAST MORE POSITIVE C IOPT (I): EARTH RADIUS OPTION C C OUTPUTS: C X,Y (R): RECTANGULAR COORDINATES IN KM C IMPLICIT NONE REAL LAT, LON, X, Y, ORGLAT, ORGLON INTEGER IOPT REAL LON1, ORGLON1, DENOM, AK REAL F, ERADEARTH, RADEARTH, ERA C DATA RADEARTH/6378.135/, F/298.26/ ! EQUITORIAL EARTH RADIUS, 1/F C ! WGS 72 MODEL VALUES LON1=MOD(LON+720.0,360.0) ORGLON1=MOD(ORGLON+720.0,360.0) C C COMPUTE LOCAL RADIUS OF THE EARTH AT CENTER OF IMAGE C ERADEARTH=6378.0 ! USE FIXED NOMINAL VALUE IF (IOPT.EQ.2) THEN ! LOCAL RADIUS ERA=(1.-1./F) ERADEARTH=RADEARTH*ERA/SQRT(ERA*ERA*COSD(ORGLAT)**2+ * SIND(ORGLAT)**2) ENDIF C DENOM=1.0+SIND(ORGLAT)*SIND(LAT)+ 1 COSD(ORGLAT)*COSD(LAT)*COSD(LON1-ORGLON1) IF (DENOM.GT.0.0) THEN AK=SQRT(2.0/DENOM) ELSE WRITE(*,*) '*** DIVISION ERROR IN LAMBERT1 ROUTINE ***' AK=1.0 ENDIF X=AK*COSD(LAT)*SIND(LON1-ORGLON1) Y=AK*(COSD(ORGLAT)*SIND(LAT)- 1 SIND(ORGLAT)*COSD(LAT)*COSD(LON1-ORGLON1)) X=X*ERADEARTH Y=Y*ERADEARTH RETURN END C C SUBROUTINE ILAMBERT1(LAT,LON,X,Y,ORGLAT,ORGLON,IOPT) C C COMPUTES THE INVERSE TRANSFORMATION FROM LAT/LON TO X/Y FOR THE C LAMBERT AZIMUTHAL EQUAL-AREA PROJECTION C C SEE "MAP PROJECTIONS USED BY THE U.S. GEOLOGICAL SURVEY" C GEOLOGICAL SURVEY BULLETIN 1532, PGS 157-173 C C FOR THIS ROUTINE, A SPHERICAL EARTH IS ASSUMED FOR THE PROJECTION. C THE ERROR WILL BE SMALL FOR SMALL-SCALE MAPS. C FOR IOPT=1 A FIXED, NOMINAL EARTH RADIUS IS USED. C FOR IOPT=2 THE LOCAL RADIUS OF THE EARTH IS USED BASED ON C THE 1972 WGS ELLIPSOID MODEL (BULLETIN PG 15). C C WRITTEN BY: DGL 7 FEB 1992 C MODIFIED: DGL 4 MAR 1995 C + INCLUDED RADIUS OF THE EARTH COMPUTED AT LOCAL REFERENCE POINT C C INPUTS: C X,Y (R): RECTANGULAR COORDINATES IN KM C ORGLAT (R): ORIGIN PARALLEL +90 TO -90 DEG WITH NORTH POSITIVE C ORGLON (R): CENTRAL MERIDIAN (LONGITUDE) 0 TO +360 DEG C OR -180 TO +180 WITH EAST MORE POSITIVE C IOPT (I): EARTH RADIUS OPTION C C OUTPUTS: C LAT (R): LATITUDE +90 TO -90 DEG WITH NORTH POSITIVE C LON (R): LONGITUDE 0 TO +360 DEG WITH EAST POSITIVE C IMPLICIT NONE REAL LAT, LON, X, Y, ORGLAT, ORGLON INTEGER IOPT REAL ORGLON1, RHO, C, T1, T2, X1, Y1 REAL F, ERADEARTH, RADEARTH, ERA C DATA RADEARTH/6378.135/, F/298.26/ ! EQUITORIAL EARTH RADIUS, 1/F C ! WGS 72 MODEL VALUES ORGLON1=MOD(ORGLON+720.0,360.0) C C COMPUTE LOCAL RADIUS OF THE EARTH AT CENTER OF IMAGE C ERADEARTH=6378.0 ! USE FIXED NOMINAL VALUE IF (IOPT.EQ.2) THEN ! LOCAL RADIUS ERA=(1.-1./F) ERADEARTH=RADEARTH*ERA/SQRT(ERA*ERA*COSD(ORGLAT)**2+ * SIND(ORGLAT)**2) ENDIF C X1=X/ERADEARTH Y1=Y/ERADEARTH RHO=X1*X1+Y1*Y1 IF (RHO.GT.0.0) THEN RHO=SQRT(RHO) C=2*ASIN(RHO*0.5) LAT=ASIND(COS(C)*SIND(ORGLAT)+Y1*SIN(C)*COSD(ORGLAT)/RHO) ELSE LAT=ORGLAT ENDIF IF (ABS(ORGLAT).NE.90.0) THEN IF (RHO.EQ.0.0) THEN LON=ORGLON1 ELSE T1=X1*SIN(C) T2=RHO*COSD(ORGLAT)*COS(C)-Y1*SIND(ORGLAT)*SIN(C) LON=ORGLON1+ATAN2D(T1,T2) ENDIF ELSE IF (ORGLAT.EQ.90.0) THEN LON=ORGLON1+ATAN2D(X1,-Y1) ELSE LON=ORGLON1+ATAN2D(X1,Y1) ENDIF LON=MOD(LON+720.0,360.0) IF (LON.GT.180.0) LON=LON-360.0 RETURN END C C SUBROUTINE POLSTER(ALON,ALAT,X1,Y1,XLAM,SLAT) C C COMPUTES THE POLAR STEROGRAPHIC TRANSFORMATION FOR A LON,LAT C INPUT OF (ALON,ALAT) WITH REFERENCE ORIGIN LON,LAT=(XLAM,SLAT). C OUTPUT IS (X1,Y1) IN KM C C ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES C AS RECEIVED FROM M. DRINKWATER (1994) C C WRITTEN BY: DGL OCT 1994 C IMPLICIT NONE REAL ALON, ALAT, X1, Y1, XLAM, SLAT C DOUBLE PRECISION E2,RE,E,T,TX,TY,CM,RHO,RLAT,SN DATA E2/0.006693883D0/ DATA RE/6378.273D0/ C E=SQRT(E2) IF (SLAT.LT.0.0) THEN SN=-1.0 RLAT=-ALAT ELSE SN=1.0 RLAT=ALAT ENDIF T=((1.D0-E*SIND(RLAT))/(1.D0+E*SIND(RLAT)))**(E*0.5D0) TY=TAND(45.0D0-0.5D0*RLAT)/T IF (SLAT.LT.0.0) THEN RLAT=-SLAT ELSE RLAT=SLAT ENDIF T=((1.D0-E*SIND(RLAT))/(1.D0+E*SIND(RLAT)))**(E*0.5D0) TX=TAND(45.0D0-0.5D0*RLAT)/T CM=COSD(RLAT)/SQRT(1.D0+E2*SIND(RLAT)**2) RHO=RE*CM*TY/TX X1= (SN*SIND(SN*ALON-XLAM))*RHO Y1=-(SN*COSD(SN*ALON-XLAM))*RHO RETURN END C C SUBROUTINE IPOLSTER(ALON,ALAT,X,Y,XLAM,SLAT) C C COMPUTES THE INVERSE POLAR STEROGRAPHIC TRANSFORMATION FOR (X,Y) C GIVEN IN KM WITH REFERENCES LON,LAT=(XLAM,SLAT). C OUTPUT LON,LAT=ALON,ALAT C C ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES C AS RECEIVED FROM M. DRINKWATER (1994). UPDATED BY D. LONG TO C IMPROVE ACCURACY USING ITERATION WITH FORWARD TRANSFORM. C C WRITTEN BY: DGL OCT 1994 C MODIFIED BY: DGL MAR 1995 C + INCORPORATED ITERATIVE UPDATE TO IMPROVE ACCURACY C IMPLICIT NONE REAL ALON, ALAT, X, Y, XLAM, SLAT C INTEGER ICNT REAL XX,YY,R,RR,RERR,A,AA,AERR,SN1 DOUBLE PRECISION E2,RE,RHO,E,E22,E23 DOUBLE PRECISION CM,CHI,T,PI2,X1,Y1,SN,SLAT1 C REAL ARCTAND EXTERNAL ARCTAND C DATA E2/0.006693883D0/ DATA RE/6378.273D0/ DATA PI2/1.570796327D0/ C C FIRST USE APPROXIMATE INVERSE CALCULATION C E=SQRT(E2) E22=E2*E2 E23=E2*E2*E2 X1=X Y1=Y RHO=X1*X1+Y1*Y1 IF (RHO.GT.0.0) RHO=SQRT(RHO) IF (RHO.LT.0.05) THEN ALON=XLAM ALAT=SIGN(90.0,SLAT) ELSE SN=1.0D0 SLAT1=SLAT IF (SLAT.LT.0.0) THEN SN=-1.0D0 SLAT1=-SLAT ENDIF CM=COSD(SLAT1)/SQRT(1.0D0-E2*SIND(SLAT1)**2) T=TAND(45.0-0.5*SLAT1)/ * ((1.D0-E*SIND(SLAT1))/(1.D0+E*SIND(SLAT1)))**(E*0.5D0) T=RHO*T/(RE*CM) CHI=PI2-2.0D0*ATAN(T) T=CHI+ * (0.5D0*E2+5.0D0*E22/24.0D0+E23/12.0D0)*SIN(2.0D0*CHI)+ * (7.0D0*E22/48.0D0+29.0D0*E23/240.0D0)*SIN(4.0D0*CHI)+ * (7.0D0*E23/120.0D0)*SIN(6.0D0*CHI) ALAT=SN*(T*90.0D0/PI2) ALON=SN*ATAN2D(SN*X1,-SN*Y1)+XLAM IF (ALON.LT.-180.0) ALON=ALON+360.0 IF (ALON.GT.180.0) ALON=ALON-360.0 ENDIF C C USING THE APPROXIMATE RESULT AS A STARTING POINT, ITERATE TO IMPROVE C THE ACCURACY OF THE INVERSE SOLUTION C SN1=1.0 IF (SLAT.LT.0.0) SN1=-1.0 A=ARCTAND(Y,X) R=SQRT(X*X+Y*Y) ICNT=0 10 CONTINUE CALL POLSTER(ALON,ALAT,XX,YY,XLAM,SLAT) RR=SQRT(XX*XX+YY*YY) RERR=SN1*(RR-R)/180.0 AA=ARCTAND(YY,XX) AERR=AA-A IF (ABS(AERR).GT.180.0) AERR=360.0-AERR C C CHECK FOR CONVERGENCE C IF ((ABS(RERR).LT.0.001.AND.ABS(AERR).LT.0.001) * .OR.ICNT.GT.9) GOTO 40 C C CONSTRAIN UPDATES C ALON=ALON+AERR IF (ABS(ALON).GT.360.0) ALON=MOD(ALON,360.0) IF (ALAT*SLAT.LT.0.0) THEN RERR=RERR*(1.0-SIND(ABS(ALAT))) IF (ABS(RERR).GT.2.0) RERR=SIGN(2.0,RERR)/FLOAT(ICNT) ENDIF ALAT=ALAT+RERR IF (ABS(ALAT).GT.90.0) ALAT=SIGN(90.0,ALAT) ICNT=ICNT+1 GOTO 10 40 CONTINUE IF (ABS(ALON).GT.360.0) ALON=MOD(ALON,360.0) RETURN END C C SUBROUTINE EASEGRID(IOPT,ALAT,ALON,THELON,THELAT,ASCALE) C C COMPUTES THE FORWARD "EASE" GRID TRANSFORM C C GIVEN A LAT,LON (ALAT,ALON) AND THE SCALE (ASCALE) THE IMAGE C TRANSFORMATION COORDINATES (THELON,THELAT) ARE COMUTED C USING THE "EASE GRID" (VERSION 1.0) TRANSFORMATION GIVEN IN FORTRAN C SOURCE CODE SUPPLIED BY NSIDC. C C THE RADIUS OF THE EARTH USED IN THIS PROJECTION IS IMBEDDED INTO C ASCALE WHILE THE PIXEL DIMENSION IN KM IS IMBEDDED IN BSCALE C THE BASE VALUES ARE: RADIUS EARTH= 6371.228 KM C PIXEL DIMEN =25.067525 KM C THEN, BSCALE = 2*(BASE_PIXEL_DIMEN) C ASCALE = 4*(BASE_RADIUS_EARTH) C C THESE VALUES ARE SET IN THE GETSIR AND PUTSIR ROUTINES C C IOPT IS EASE TYPE: IOPT=11=NORTH, IOPT=12=SOUTH, IOPT=13=CYLINDRICAL C C WRITTEN BY: DGL 4 MAR 1995 C IMPLICIT NONE INTEGER IOPT REAL ALAT, ALON, THELON, THELAT, ASCALE REAL PI2 DATA PI2/1.57079633/ ! Pi/2 C IF (IOPT.EQ.11) THEN ! EASE GRID NORTH THELON= ASCALE*SIND(ALON)*SIND(45.0-0.5*ALAT) THELAT=-ASCALE*COSD(ALON)*SIND(45.0-0.5*ALAT) ELSE IF (IOPT.EQ.12) THEN ! EASE GRID SOUTH THELON=ASCALE*SIND(ALON)*COSD(45.0-0.5*ALAT) THELAT=ASCALE*COSD(ALON)*COSD(45.0-0.5*ALAT) ELSE IF (IOPT.EQ.13) THEN ! EASE CYLINDRICAL THELON=ASCALE*PI2*ALON*COSD(30.0)/90.0 THELAT=ASCALE*SIND(ALAT)*COSD(30.0) ENDIF RETURN END C C SUBROUTINE IEASEGRID(IOPT,ALON,ALAT,THELON,THELAT,ASCALE) C C COMPUTES THE INVERSE "EASE" GRID TRANSFORM C C GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND C THE SCALE (ASCALE) THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED C USING THE "EASE GRID" (VERSION 1.0) TRANSFORMATION GIVEN IN FORTRAN C SOURCE CODE SUPPLIED BY NSIDC C C IOPT IS EASE TYPE: IOPT=11=NORTH, IOPT=12=SOUTH, IOPT=13=CYLINDRICAL C C THE RADIUS OF THE EARTH USED IN THIS PROJECTION IS IMBEDDED INTO C ASCALE WHILE THE PIXEL DIMENSION IN KM IS IMBEDDED IN BSCALE C THE BASE VALUES ARE: RADIUS EARTH= 6371.228 KM C PIXEL DIMEN =25.067525 KM C THEN, BSCALE = 2*(BASE_PIXEL_DIMEN) C ASCALE = 4*(BASE_RADIUS_EARTH) C C THESE VALUES ARE SET IN THE GETSIR AND PUTSIRROUTINES C C THE INVERSE ALGORITHM IS SPECIFIED BY THE EASE GRID DEFINITION. C IT IS ACTUALLY ONLY AN APPROXIMATION TO THE MODIFIED POLAR C STEREOGRAPHIC PROJECTION USED IN THE EASE GRID DEFINITION. C C CODE WRITTEN BY: DGL 4 MAR 1995 C IMPLICIT NONE INTEGER IOPT REAL ALON, ALAT, THELON, THELAT, ASCALE REAL PI2, X1, Y1,TEMP DATA PI2/1.57079633/ ! Pi/2 REAL ARCTAND ! EXTERNAL FUNCTION C X1=THELON Y1=THELAT IF (IOPT.EQ.11) THEN ! EASE GRID NORTH ALON=ARCTAND(X1,-Y1) IF (ABS(SIND(ALON)).GT.ABS(COSD(ALON))) THEN TEMP=(X1/SIND(ALON))/ASCALE ELSE TEMP=(-Y1/COSD(ALON))/ASCALE ENDIF IF (ABS(TEMP).LE.1.0) THEN ALAT=90.0-2.*ASIND(TEMP) ELSE ALAT=-90.0+2.*ASIND(1./TEMP) WRITE (*,*) '*** ERROR IN EASE INVERSE SINE ***',TEMP C ALAT=SIGN(90.0,TEMP) ENDIF ELSE IF (IOPT.EQ.12) THEN ! EASE GRID SOUTH ALON=ARCTAND(X1,Y1) IF (ABS(COSD(ALON)).GT.ABS(SIND(ALON))) THEN TEMP=(X1/COSD(ALON))/ASCALE ELSE TEMP=(Y1/SIND(ALON))/ASCALE ENDIF IF (ABS(TEMP).LE.1.0) THEN ALAT=90.0-2.*ACOSD(TEMP) ELSE WRITE (*,*) '*** ERROR IN EASE INVERSE COSINE ***',TEMP ALAT=0.0 ENDIF ELSE IF (IOPT.EQ.13) THEN ! EASE CYLINDRICAL ALON=((X1/ASCALE)/COSD(30.0))*90.0/PI2 TEMP=(Y1/COSD(30.0))/ASCALE IF (ABS(TEMP).LE.1.0) THEN ALAT=ASIND(TEMP) ELSE WRITE (*,*) '*** ERROR IN EASE INVERSE SINE ***',TEMP ALAT=SIGN(90.0,TEMP) ENDIF ENDIF RETURN END C C REAL FUNCTION ARCTAND(Y,X) IF (Y.EQ.0.0.AND.X.EQ.0.0) THEN ARCTAND=0.0 ELSE ARCTAND=ATAN2D(Y,X) ENDIF RETURN END