PROGRAM SIR_STATS C C THIS PROGRAM SUMMARIZES THE PIXEL VALUE STATISTICS C OF A SUB REGION FROM A SIR FILE C C WRITTEN BY: DGL JUL 2005 + based on sir_extractregion.f C C Uses command line arguments to a fortran program via IARGC and GETARG C Not all compilers/operating systems support this. If not supported C Comment out the CALL GETARG lines and the NARG=IARGC() line C C READVAL1 and IREADVAL1 are I/O routines to simplify input. C They can be replaced with READ (*,*) statements C C LINK to SIR library routines C CHARACTER*180 FNAME CHARACTER*190 NAME,LINE CHARACTER TITLE2*80 LOGICAL LATLON C C IMAGE DISPLAY VARIABLES C PARAMETER (MAXSIZE=16000000) REAL STVAL(MAXSIZE) REAL STVAL2(MAXSIZE) 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 BEGIN PROGRAM C WRITE (*,5) 5 FORMAT(/'SIR_STATS -- summarize statistics from subregion in SIR file '/) NARG=0 C C compiler dependent number of arguments routine NARG=IARGC() ! GET NUMBER OF COMMAND LINE ARGUMENTS C IF (NARG.LT.1) THEN WRITE (*,*) 'Command line: sirextregion in_file P/L LLx LLy URx URy out_file' WRITE (*,*) ' in_file: input sir file' WRITE (*,*) ' P/L: pixel (P) or lat/lon (L) coordinates' WRITE (*,*) ' LLx,LLy: lower-left corner (L=lon,lat)' WRITE (*,*) ' URx,URy: upper-right corner (L=lon,lat)' WRITE (*,*) ' out_file: output sir file' ENDIF C IF (NARG.GT.0) THEN C compiler dependent routine to fill FNAME with first run argument CALL GETARG(1,FNAME) ELSE WRITE (*,*) 'Source file to read from?' READ(*,15) FNAME 15 FORMAT(A80) ENDIF C C READ INPUT IMAGE FILE C WRITE(*,16) FNAME(1:LENGTH(FNAME)) 16 FORMAT(' Reading file "',A,'"') MAXDES=LEN(DESCRIP) IU=1 CALL READSIRHEAD3(FNAME,IU,IERR,NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IXDEG_OFF,IYDEG_OFF,IDEG_SC,ISCALE_SC,IA0_OFF,IB0_OFF,IO_SC, * 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 ECHO INPUT FILE HEADER TO USER C WRITE (*,*) 'Input file header information:' CALL PRINTHEAD3(NHEAD,NDES,NHTYPE,IDATATYPE, * NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE, * IXDEG_OFF,IYDEG_OFF,IDEG_SC,ISCALE_SC,IA0_OFF,IB0_OFF,IO_SC, * 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 OUTPUT HEADER C IF (NHTYPE.EQ.1) THEN ! OLD STYLE HEADER 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 SMIN=MIN(SMIN,STVAL(INDEX)) SMAX=MAX(SMAX,STVAL(INDEX)) END DO END DO ANODATA=SMIN VMIN=SMIN VMAX=SMAX IF (ABS(SMIN+32.0).LT.2) THEN ANODATA=-32.0 VMIN=-32.0 ELSE IF (ABS(SMIN+3.2).LT.2) THEN ANODATA=-3.2 VMIN=-3.0 ENDIF IF (ABS(SMAX).LT.2) THEN VMAX=0.0 ENDIF ENDIF C LATLON=.TRUE. IF (NARG.GE.2) THEN CALL GETARG(2,LINE) IF (LINE(1:1).EQ.'P'.OR.LINE(1:1).EQ.'p') LATLON=.FALSE. ELSE WRITE (*,514) 514 FORMAT('Pixel (P) or Lat/Lon (L [def]) coordinates: ',$) READ(*,15) LINE IF (LINE(1:1).EQ.'P'.OR.LINE(1:1).EQ.'p') LATLON=.FALSE. ENDIF C IF (LATLON) THEN C IF (NARG.GE.3) THEN CALL GETARG(3,LINE) READ(LINE,'(F)') ALON1 ELSE WRITE (*,516) 516 FORMAT('Lower-left corner longitude (X) ',$) C following can be replaced with READ(*,*) ALON1 ALON1=READVAL1(ALON1,IEND) IF (IEND.LT.0) STOP ENDIF C IF (NARG.GE.4) THEN CALL GETARG(4,LINE) READ(LINE,'(F)') ALAT1 ELSE WRITE (*,515) 515 FORMAT('Lower-left corner latitude (Y) ',$) C following can be replaced with READ(*,*) ALA1 ALAT1=READVAL1(ALAT1,IEND) IF (IEND.LT.0) STOP ENDIF C IF (NARG.GE.5) THEN CALL GETARG(5,LINE) READ(LINE,'(F)') ALON2 ELSE WRITE (*,518) 518 FORMAT('Upper-right corner longitude (X) ',$) C following can be replaced with READ(*,*) ALON2 ALON2=READVAL1(ALON2,IEND) IF (IEND.LT.0) STOP ENDIF C IF (NARG.GE.6) THEN CALL GETARG(6,LINE) READ(LINE,'(F)') ALAT2 ELSE WRITE (*,517) 517 FORMAT('Upper-right corner latitude (Y) ',$) C following can be replaced with READ(*,*) ALAT2 ALAT2=READVAL1(ALAT2,IEND) IF (IEND.LT.0) STOP ENDIF C C CONVERT LAT, LON CORNERS TO IMAGE PIXELS C CALL LATLON2PIX(ALON1,ALAT1,X1,Y1,THELON,THELAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) CALL F2IPIX(X1,Y1,IX1,IY1,NSX,NSY) WRITE (*,*) 'Input lon,lat: ',ALON1,ALAT1,' -> ',IX1,IY1 IF (IX1.EQ.0.OR.IY1.EQ.0) THEN WRITE (*,*) '** POINT OUTSIDE OF IMAGE ***' IX1=MIN(1,IX1) IY1=MIN(1,IY1) ENDIF C CALL LATLON2PIX(ALON2,ALAT2,X2,Y2,THELON,THELAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) CALL F2IPIX(X2,Y2,IX2,IY2,NSX,NSY) WRITE (*,*) 'Input lon,lat: ',ALON2,ALAT2,' -> ',IX2,IY2 IF (IX2.EQ.0.OR.IY2.EQ.0) THEN WRITE (*,*) '** POINT OUTSIDE OF IMAGE ***' IX2=MIN(1,IX2) IY2=MIN(1,IY2) ENDIF C ELSE C IF (NARG.GE.3) THEN CALL GETARG(3,LINE) READ(LINE,'(I)') IX1 ELSE WRITE (*,1516) 1516 FORMAT('Lower-left corner pixel (X) ',$) C following can be replaced with READ(*,*) IX1 IX1=IREADVAL1(1,IEND) IF (IEND.LT.0) STOP ENDIF C IF (NARG.GE.4) THEN CALL GETARG(4,LINE) READ(LINE,'(I)') IY1 ELSE WRITE (*,1515) 1515 FORMAT('Lower-left corner pixel (Y) ',$) C following can be replaced with READ(*,*) IY1 IY1=IREADVAL1(1,IEND) IF (IEND.LT.0) STOP ENDIF C IF (NARG.GE.5) THEN CALL GETARG(5,LINE) READ(LINE,'(I)') IX2 ELSE WRITE (*,1518) 1518 FORMAT('Upper-right corner pixel (X) ',$) C following can be replaced with READ(*,*) IX2 IX2=IREADVAL1(NSX,IEND) IF (IEND.LT.0) STOP ENDIF C IF (NARG.GE.6) THEN CALL GETARG(6,LINE) READ(LINE,'(I)') IY2 ELSE WRITE (*,1517) 1517 FORMAT('Upper-right corner pixel (Y) ',$) C following can be replaced with READ(*,*) IY2 IY2=IREADVAL1(NSY,IEND) IF (IEND.LT.0) STOP ENDIF ENDIF C C CLIP SUBREGION TO BE WITHIN IMAGE AND SORT CORNERS C IF (IX1.LT.1) IX1=1 IF (IY1.LT.1) IY1=1 IF (IX2.LT.1) IX2=1 IF (IY2.LT.1) IY2=1 C IF (IX1.GT.NSX) IX1=NSX IF (IY1.GT.NSY) IY1=NSY IF (IX2.GT.NSX) IX2=NSX IF (IY2.GT.NSY) IY2=NSY C IX=MIN(IX1,IX2) IX2=MAX(IX1,IX2) IX1=IX IY=MIN(IY1,IY2) IY2=MAX(IY1,IY2) IY1=IY NSX2=IX2-IX1+1 NSY2=IY2-IY1+1 C X=IX1 Y=IY1 CALL PIX2LATLON(X,Y,THELON,THELAT,ALON,ALAT, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) C X=IX2 Y=IY2 CALL PIX2LATLON(X,Y,THELON,THELAT,ALON2,ALAT2, $ IOPT,XDEG,YDEG,ASCALE,BSCALE,A0,B0) C WRITE (*,*) 'Image size: (in) ',NSX,NSY,' (out) ',NSX2,NSY2,' ',NSX2*NSY2 WRITE (*,*) 'Lower-Left pixel location: ',IX1,IY1,' Lon, Lat: ',ALON,ALAT WRITE (*,*) 'Upper-Right pixel location:',IX2,IY2,' Lon, Lat: ',ALON2,ALAT2 C IF (NSX2.LT.1.OR.NSY2.LT.1) THEN WRITE (*,*) ' *** Error Extracting SubRegion: ZERO SIZE ',NSX2,NSY2 STOP ENDIF C C CALCULATE IMAGE STATISTICS C NCNT=0 AVE=0.0 RMS=0.0 DMAX1=-1.e25 DMIN1=1.e25 DO IY=1,NSY2 DO IX=1,NSX2 I1=(IY-1+IY1-1)*NSX+(IX+IX1-1) IF (STVAL(I1).GT.ANODATA+0.0001) THEN NCNT=NCNT+1 AVE=STVAL(I1)/FLOAT(NCNT)+AVE*FLOAT(NCNT-1)/FLOAT(NCNT) RMS=STVAL(I1)*STVAL(I1)/FLOAT(NCNT)+RMS*FLOAT(NCNT-1)/FLOAT(NCNT) IF (STVAL(I1).GT.DMAX1) DMAX1=STVAL(I1) IF (STVAL(I1).LT.DMIN1) DMIN1=STVAL(I1) ENDIF END DO END DO C C SUMMARIZE STATS C IF (NCNT.EQ.0) THEN WRITE (*,*) ' *** No valid pixels in area' STOP ENDIF IF (RMS.GT.0.0) THEN STD=RMS-AVE*AVE IF (STD.GT.0.0) STD=SQRT(STD) RMS=SQRT(RMS) ENDIF C WRITE (*,*) ' pix cnt: ',NCNT,' of ',NSX2*NSY2 WRITE (*,*) ' Ave, std: ',AVE,STD WRITE (*,*) ' Min, Max: ',DMIN1,DMAX1 STOP END