	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
