PROGRAM SIR2BYTE C C THIS PROGRAM ILLUSTRATES THE READING OF A FILE IN THE BYU MERS SIR C FILE FORMAT. IT READS THE FILE AND C WRITES OUT A BYTE-ARRAY COPY OF THE IMAGE C C WRITTEN BY D.LONG: MARCH 1997 C CHARACTER*70 FNAME,ONAME C C IMAGE STORAGE VARIABLES C PARAMETER (MAXSIZE=4500000) ! MAXIMUM PIXELS REAL STVAL(MAXSIZE) ! READ ARRAY BYTE B(MAXSIZE) ! OUTPUT BYTE ARRAY 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 READ IMAGE CONTENTS C IF (IERR.GE.0) THEN CALL READSIRF(IU,IERR,NHEAD,NHTYPE,IDATATYPE,STVAL,NSX,NSY, * IOFF,ISCALE,SMIN,SMAX,NCNT,ANODATA,VMIN,VMAX) ELSE WRITE (*,*) '*** ERROR OPENING FILE',IERR ENDIF IF (IERR.LT.0) STOP 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 MAX/MIN OF DATA ARRAY C SMIN=1.E25 SMAX=-1.E25 DO NN1=0,NSY-1 DO NN2=1,NSX INDEX = NN1*NSX + NN2 IF (STVAL(INDEX).NE.ANODATA) THEN SMIN=MIN(SMIN,STVAL(INDEX)) SMAX=MAX(SMAX,STVAL(INDEX)) ENDIF END DO END DO C SMIN1=SMIN SMAX1=SMAX C WRITE (*,*) WRITE (*,*) 'Array min, max values:',SMIN1,SMAX1 C WRITE (*,34) 'Saturation Max?' 34 FORMAT(X,'New ',A,X,$) READ (*,*) SMAX WRITE (*,34) 'Saturation Min?' READ (*,*) SMIN C C SCALE ARRAY AND PUT INTO BYTE ARRAY C C NOTE THAT SIR IMAGE DATA HAS THE ORIGIN AT THE LOWER-LEFT OF THE C IMAGE. THE BYTE ARRAY HAS THE ORIGIN AT THE UPPER-LEFT. C SS=(SMAX-SMIN)/255.0 IF (SS.GT.0.0) THEN SS=1.0/SS ELSE SS=1.0 ENDIF DO NN1=0,NSY-1 DO NN2=1,NSX INDEX = NN1*NSX + NN2 INDEX2= (NSY-NN1)*NSX + NN2 VAL=STVAL(INDEX) IVAL=(VAL-SMIN)*SS IF (IVAL.LT.0) IVAL=0 IF (IVAL.GT.255) IVAL=255 IF (IVAL.GT.127) IVAL=IVAL-256 B(INDEX2)=IVAL END DO END DO C C GET OUTPUT FILE NAME FOR BYTE FILE C WRITE(*,25) 25 FORMAT(' Enter BYTE Output File Name: ',$) READ(*,20) ONAME WRITE(*,*) 'BYTE Out="',ONAME(1:LENGTH(ONAME)),'"' C C NOW WRITE OUT BYTE FILE C WRITE(*,*) 'Output Image Size: (X x Y)',NSX,NSY NS=NSX*NSY C OPEN(UNIT=1,FILE=ONAME,STATUS='NEW',FORM='UNFORMATTED') C WRITE (1) (B(I),I=1,NS) CLOSE (1) 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)',ascale,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 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 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 ***************************************************************************