/******************************************************************************
* 
* grid2grid 
*
*   Regrids the input grid to specified output grid resolution, averaging
*   or summing the input grid values (OR, generating a grid of maximum or 
*   minimum values per bin).  Weighted averaging is used to handle data from 
*   input grid cells that overlap more than one output grid cell. 
*
* LAST UPDATE: 21 May 1999 
*
* USAGE:
*   "grid2grid" prompts the user for input  
*
* ARGUMENTS:
*    None.
* 
* ORIGINAL AUTHOR:
*   Wes Berg, CIRES/CDC, 1995(?) 
*
* NOTES:
*   (1) Use -lm flag to compile (e.g., cc grid2grid.c -lm). 
*   (2) Currently hard-coded to a limit of 1000 x 1000 output grid size. 
*   (3) Programmer's note: a single grid array is defined to store output 
*       values for any of the 4 output options (mean, sum, max, min). 
*       Defining separate grids requires allocating more memory than is
*       normally available (i.e., it will result in segmentation faults).
*
* MODIFICATIONS:
*   Don Anderson  05/20/98  Modified to work on any input data type 
*   Don Anderson  05/20/98  Modified to output regridded means OR sums  
*   Don Anderson  05/11/99  Corrected to create non-rounded float output  
*   Don Anderson  05/21/99  Added option to output regridded max or min values  
*
******************************************************************************/
#include <stdio.h>
#include <math.h>
#define PMODE 0644                /* Output file permissions             */
#define SMAX 1000                 /* Maximum number of samples in line   */ 
#define CUTOFF 0.5                /* Added to floats to calc nearest int */ 
#define NLMAX 1000                /* Max lines in output grid            */ 
#define NSMAX 1000                /* Max samples in output grid          */ 
#define NBMAX 5000                /* Max number of bands                 */ 

void main()
 {
 float  clat, clon ;               /* Output pixel center lat, lon         */
 int    cut ;                      /* Flag: cut off values at 0? (1=y 0=n) */ 
 unsigned char  *databuff1 ;       /* Byte data buffer                     */ 
 short *databuff2 ;                /* Short data buffer                    */ 
 long  *databuff3 ;                /* Long data buffer                     */ 
 float *databuff4 ;                /* Float data buffer                    */ 
 int    datasize ;                 /* Bytes per input data value           */
 int    datatype ;                 /* Indata type 1=byt 2=sht 3=lng 4=flt  */ 
 int    first ;                    /* Flag for first interpolation loop    */ 
 static float grid[NLMAX][NSMAX] ; /* Grid of summed weighted data values  */ 
 int    i ;                        /* Line counter                         */
 int    ilat[4] ;                  /* Output xdim cell indices (integers)  */
 int    ilon[4] ;                  /* Output ydim cell indices (integers)  */ 
 int    inf ;                      /* Input file handle                    */ 
 char   infile[80] ;               /* Input file string                    */ 
 short  ival ;                     /* Output value rounded to nearest int  */ 
 int    j ;                        /* Sample counter                       */ 
 int    k ;                        /* Generic counter(?)                   */
 int    l ;                        /* Counter for 4 output cell corners    */ 
 float  lat[4] ;                   /* Input cell four corner latitudes     */ 
 float  lon[4] ;                   /* Input cell four corner longitudes    */ 
 int    m ;                        /* Generic counter                      */
 int    mflag ;                    /* Do-not-interpolate if flag == 1      */ 
 float  min[NBMAX], max[NBMAX] ;   /* Min, max input values by band        */ 
 float  min2[NBMAX], max2[NBMAX] ; /* Min, max output values by band       */ 
 float  minlat1, maxlat1 ;         /* Min/max latitude, input file         */ 
 float  minlon1, maxlon1 ;         /* Min/max longitude, input file        */ 
 float  minlat2, maxlat2 ;         /* Min/max latitude, output file        */
 float  minlon2, maxlon2 ;         /* Min/max longitude, output file       */ 
 double atof() ;                   /* Function to convert string to float  */
 int    n ;                        /* Number of values to not interpolate  */
 int    nb ;                       /* Band counter                         */ 
 int    nbnd1 ;                    /* Number of input file bands           */ 
 int    nl ;                       /* Number of input lines                */ 
 int    ns ;                       /* Number of input samples              */ 
 int    nlout ;                    /* Number of output lines               */
 int    nsout ;                    /* Number of output samples             */
 int    nsp[20] ;                  /* Array of values to not interpolate   */ 
 float  nweights[NLMAX][NSMAX] ;   /* No-interp summed weighting factors   */
 int    opf ;                      /* Output file handle                   */ 
 char   outfile[80] ;              /* Output file string                   */
 unsigned char  *outbuf1 ;         /* Byte data buffer                     */ 
 short *outbuf2 ;                  /* Short data buffer                    */ 
 long  *outbuf3 ;                  /* Long data buffer                     */ 
 float *outbuf4 ;                  /* Float data buffer                    */ 
 int    outflag ;                  /* Output: 0=sum, 1=mean, 2=max, 3=min  */ 
 float  rlat ;                     /* Input cell center latitude           */ 
 float  rlon ;                     /* Input cell center longitude          */ 
 char   s[10] ;                    /* Generic string holder                */ 
 float  scale ;                    /* Input data scale factor              */ 
 int    sflag ;                    /* 1=matching ilat/ilon pair            */ 
 float  sp ;                       /* Non-interpolated pixel value?        */
 float  val ;                      /* Pixel value                          */
 float  weights[NLMAX][NSMAX] ;    /* Summed weighting factors             */
 float  wt ;                       /* Infile cell weighting factor (max 1) */ 
 float  wtx ;                      /* X weighting for averaging (max is 1) */ 
 float  wty ;                      /* Y weighting for averaging (max is 1) */ 
 float  wtmax ;                    /* Maximum weight                       */ 
 float  xdist ;                    /* Critical latitude distance (degs) from 
                                      output to input cell centers, beyond 
                                      which in-cell overlaps next out-cell */  
 float  x1size ;                   /* Input x cell dimension in degrees    */
 float  x2size ;                   /* Output x cell dimension in degrees   */
 float  ydist ;                    /* Critical longitude distance (degs) from 
                                      output to input cell centers, beyond 
                                      which in-cell overlaps next out-cell */  
 int    yval ;                     /* Non-scaled pixel value as int        */
 float  y1size ;                   /* Input y cell dimension in degrees    */
 float  y2size ;                   /* Output y cell dimension in degrees   */

 wtmax = 0.0 ; 

/***************************/
/*  open input data file   */
/***************************/

 printf (" Please enter the name of the input file: ");
 gets (infile);
 while ( (inf = open (infile,0) ) == -1) 
   {
   printf(" Unable to open input file, please re-enter! \n");
   printf(" Please enter the name of the input file: ");
   gets(infile);
   }
 printf (" Please enter the number of lines in the input file: ");
 gets (s) ;
 nl = atoi (s) ;
 printf (" Please enter the number of samples in the input file: ");
 gets (s);
 ns = atoi (s);

 printf (" Please enter the input data type (b=byte, s=short, l=long, f=float): ");
 gets (s);
 if ( ! strcmp ( s, "b") || ! strcmp ( s, "B" ) )
   {
   datatype = 1 ;
   datasize = 1 ;
   databuff1 = ( unsigned char * ) malloc (  SMAX * sizeof ( char ) ) ;
   outbuf1 = ( unsigned char * ) malloc ( SMAX * sizeof ( char ) ) ;
   }
 else if ( ! strcmp ( s, "s") || ! strcmp ( s, "S" ) )
   {
   datatype = 2 ;
   datasize = 2 ;
   databuff2 = ( short * ) malloc ( SMAX * sizeof ( short ) ) ;
   outbuf2 = ( short * ) malloc ( SMAX * sizeof ( short ) ) ;
   }
 else if ( ! strcmp ( s, "l") || ! strcmp ( s, "L" ) )
   {
   datatype = 3 ;
   datasize = 4 ;
   databuff3 = ( long * ) malloc ( SMAX * sizeof ( long ) ) ;
   outbuf3 = ( long * ) malloc ( SMAX * sizeof ( long ) ) ;
   }
 else if ( ! strcmp ( s, "f") || ! strcmp ( s, "F" ) )
   {
   datatype = 4 ;
   datasize = 4 ;
   databuff4 = ( float * ) malloc ( SMAX * sizeof ( float ) ) ;
   outbuf4 = ( float * ) malloc ( SMAX * sizeof ( float ) ) ;
   }
 else 
   {
   printf (" Invalid data type option! ") ;
   exit (1) ;
   } 

 if ( ( ! databuff1 || ! outbuf1 ) &&  
      ( ! databuff2 || ! outbuf2 ) &&
      ( ! databuff3 || ! outbuf3 ) &&
      ( ! databuff4 || ! outbuf4 ) )  
   {
   printf (" Error allocating memory! ") ;
   exit (1) ;
   } 

 printf (" Please enter the number of bands in the input file: ");
 gets (s);
 nbnd1 = atoi (s);
 printf (" Please enter the scale factor for the input data file: ");
 gets (s);
 scale = atof (s);

/***************************/
/*  open output data file  */
/***************************/

 printf (" Please enter the name of the output file: ");
 gets (outfile);
 while ( ( opf = creat ( outfile, PMODE ) ) == -1 ) 
   {
   printf(" Unable to open output file, please re-enter! \n");
   printf(" Please enter the name of the output file: ");
   gets ( outfile ) ;
   }
 printf ( " Please enter the number of lines for the output image: ");
 gets (s);
 nlout = atoi (s);
 printf(" Please enter the number of samples for the output image: ");
 gets (s);
 nsout = atoi(s);
 printf (" Please enter the number of values to not interpolate: ");
 gets (s);
 n = atoi (s);
 for ( i=0; i<n; i++ ) 
   {
   printf (" Please enter the value to remain constant (integer): ");
   gets (s);
   nsp[i] = atoi (s);
   }
 sp = -5.0 ;

 printf (" Does the input file region differ from 360 degrees longitude, 60N/60S (y/n)? : ");
 gets (s);
 if ( ! strcmp ( s, "Y") || ! strcmp ( s, "y" ) ) 
   {
   printf (" Please enter the minimum latitude for the input file: ");
   gets (s);
   minlat1 = atof (s);
   printf (" Please enter the maximum latitude for the input file: ");
   gets (s);
   maxlat1 = atof(s);
   printf (" Please enter the minimum longitude for the input file: ");
   gets (s);
   minlon1 = atof (s);
   printf (" Please enter the maximum longitude for the input file: ");
   gets (s);
   maxlon1 = atof (s);
   }
 else
   {
   minlat1 = -60.0 ; 
   maxlat1 = 60.0 ; 
   minlon1 = -180.0 ; 
   maxlon1 = 180.0 ; 
   }

 printf (" Do you wish to define an output file region that is different from the input (y/n)? : ");
 gets (s);
 if ( ! strcmp ( s, "Y") || ! strcmp ( s, "y" ) ) 
   {
   printf (" Please enter the minimum latitude for the output file: ");
   gets(s);
   minlat2 = atof (s);
   printf (" Please enter the maximum latitude for the output file: ");
   gets (s);
   maxlat2 = atof (s);
   printf (" Please enter the minimum longitude for the output file: ");
   gets (s);
   minlon2 = atof (s);
   printf (" Please enter the maximum longitude for the output file: ");
   gets (s);
   maxlon2 = atof (s);
   }
 else
   {
   minlat2 = minlat1 ; 
   maxlat2 = maxlat1 ; 
   minlon2 = minlon1 ; 
   maxlon2 = maxlon1 ; 
   }

 printf (" Do you want output to be regridded sums, averages, maxima or minima (s/a/x/n)? ") ;
 gets (s) ;
 if ( s[0] == 'S' || s[0] == 's' ) outflag = 0 ;
 else if ( s[0] == 'M' || s[0] == 'm' ) outflag = 1 ;
 else if ( s[0] == 'X' || s[0] == 'x' ) outflag = 2 ;
 else if ( s[0] == 'N' || s[0] == 'n' ) outflag = 3 ;
 else outflag = 1 ; 

 printf (" Do you wish to cut off values at 0.0 (y/n)? ");
 gets (s);
 if ( s[0] == 'y' || s[0] == 'Y' )
   cut = 1;
 else
   cut = 0;

 x1size = ( maxlon1 - minlon1 ) / (float) (ns);
 y1size = ( maxlat1 - minlat1 ) / (float) (nl);
 x2size = ( maxlon2 - minlon2 ) / (float) (nsout);
 y2size = ( maxlat2 - minlat2 ) / (float) (nlout);
 xdist = ( x1size / 2.0 + x2size / 2.0 );
 ydist = ( x1size / 2.0 + x2size / 2.0 );

/*********************************/
/*  read in the input data file  */
/*********************************/

/** Loop thru input bands **/ 

 for ( nb=0; nb<nbnd1; nb++ ) 

   {

/** Initialize grid, weights and nweights arrays to zeros **/

   for ( i=0; i<nlout; i++ ) 
     {
     for ( j=0; j<nsout; j++ ) 
       {
       if ( outflag == 0 || outflag == 1 ) grid[i][j] = 0.0 ;
       else if ( outflag == 2 ) grid[i][j] = -99999.0 ;
       else if ( outflag == 3 ) grid[i][j] = 99999.0 ;
       else grid[i][j] = 0.0 ;
       weights[i][j] = 0.0 ;
       nweights[i][j] = 0.0 ;
       }
     }

/** Loop thru input lines **/

   first = 0 ;
   for ( i=0; i<nl; i++ ) 
     {
     rlat = minlat1 + ( (float) i + 0.5 ) * y1size ;
     if      ( datatype == 1 ) read ( inf, databuff1, (size_t) ns * datasize ) ;
     else if ( datatype == 2 ) read ( inf, databuff2, (size_t) ns * datasize ) ;
     else if ( datatype == 3 ) read ( inf, databuff3, (size_t) ns * datasize ) ;
     else if ( datatype == 4 ) read ( inf, databuff4, (size_t) ns * datasize ) ;

/** Loop thru input samples **/

     for ( j=0; j<ns; j++ ) 
       {
       rlon = minlon1 + ( (float) j + 0.5 ) * x1size ;

/** Check to see if input cell edge falls w/in output region **/ 

       if ( rlat + y1size / 2.0 > minlat2 &&
            rlat - y1size / 2.0 < maxlat2 &&
            rlon + x1size / 2.0 > minlon2 &&
            rlon - x1size / 2.0 < maxlon2 ) 
         {
         mflag = 0 ;

/** If so, extract the scaled (val) and original (yval) data values **/
 
         if      ( datatype == 1 ) val = scale * (float) databuff1[j] ;
         else if ( datatype == 2 ) val = scale * (float) databuff2[j] ;
         else if ( datatype == 3 ) val = scale * (float) databuff3[j] ;
         else if ( datatype == 4 ) val = scale * (float) databuff4[j] ;

         yval = ( int ) ( val / scale ) ; 

/** Compare against list of values for which no interpolation wanted **/  
/** Set val to -5 if not to be interpolated                          **/

         for ( k=0; k<n; k++ ) 
           {
           if ( yval == nsp[k] )
             {
             val = sp ;
             mflag = 1 ;
             }
           }

/** Keep track of the minimum and maximum values by band **/   

         if ( ! mflag ) 
           {
           if ( first == 0 ) 
             {
             min[nb] = val ;
             max[nb] = val ;
             first = 1;
             }
           else 
             {
             if ( val > max[nb] ) max[nb] = val ;
             if ( val < min[nb] ) min[nb] = val ;
             }
           }

/** Calculate latitudes & longitudes of the four corners of the input cell **/
/** Indices: 0=upper right, 1=upper left, 2=lower right, 3=lower left      **/

         lat[0] = rlat + y1size / 2.0 ;
         lon[0] = rlon + x1size / 2.0 ;
         lat[1] = rlat + y1size / 2.0 ;
         lon[1] = rlon - x1size / 2.0 ;
         lat[2] = rlat - y1size / 2.0 ;
         lon[2] = rlon + x1size / 2.0 ;
         lat[3] = rlat - y1size / 2.0 ;
         lon[3] = rlon - x1size / 2.0 ;


/** Calculate the equivalent output cell x,y index location for each corner **/

         for ( l=0; l<4; l++ ) 
           {
           ilat[l] = (int) ( ( lat[l] - minlat2 ) / y2size );
           if ( ilat[l] < 0 ) ilat[l] = 0;
           if ( ilat[l] > nlout-1 ) ilat[l] = nlout-1 ;
           ilon[l] = (int) ( ( lon[l] - minlon2 ) / x2size ) ;
           if ( ilon[l] < 0 ) ilon[l] = 0;
           if ( ilon[l] > nsout-1 ) ilon[l] = nsout-1 ;


/** Check for matching output cell x,y locations for different corners **/

           sflag = 0 ;
           for ( m=0; m<l; m++ ) 
             {
             if ( ilat[l] == ilat[m] && ilon[l] == ilon[m] ) 
               {
               sflag = 1 ;
               break ;
               }
             }

/** Determine weighting factors if not already calced for this ilat/ilon **/

           if ( ! sflag ) 
             {
             clat = minlat2 + ( (float) ilat[l] + 0.5 ) * y2size ;
             clon = minlon2 + ( (float) ilon[l] + 0.5 ) * x2size ;
             wtx = ( xdist - fabs ( clon - rlon ) ) / x1size ;
             if ( wtx > 1.0 ) wtx = 1.0 ;
             wty = ( ydist - fabs ( clat - rlat ) ) / y1size ;
             if ( wty > 1.0 ) wty = 1.0 ;
             wt = wtx * wty ;

/** Add the weighting factor to the weights matrix **/ 

             weights[ilat[l]][ilon[l]] += wt ;

/** Sum the weighted values (val*wt) values in grid matrix  **/ 
/** (unless no-interp flag, then sum wt in nweights matrix) **/
/** Also, record new max and min values in grid matrices    **/ 

             if ( ! mflag )
               {
               if ( outflag == 0 || outflag == 1 ) 
                 grid[ilat[l]][ilon[l]] += val * wt ;
               if ( outflag == 2 )
                 if ( val > grid[ilat[l]][ilon[l]] ) 
                   grid[ilat[l]][ilon[l]] = val ;
               if ( outflag == 3 )
                 if ( val < grid[ilat[l]][ilon[l]] ) 
                   grid[ilat[l]][ilon[l]] = val ;
               }
             else
               nweights[ilat[l]][ilon[l]] += wt ;
             if ( weights[ilat[l]][ilon[l]] > wtmax ) 
               wtmax = weights[ilat[l]][ilon[l]] ;
             }
           }
         }
       }
     }

/*********************************/
/*  Perform interpolation        */
/*********************************/

   first = 0 ;

   for ( i=0; i<nlout; i++ ) 
     {
     for ( j=0; j<nsout; j++ ) 
       {

       /** If all weightings from no-interp values, keep same output value **/

       if ( nweights[i][j] / weights[i][j] == 1.0 ) 
         {
         val = sp ;
         ival = nsp[0] ;
         }
        
       /** Otherwise, evaluate actual data values **/

       else 

         {

         if ( outflag == 0 ) val = grid[i][j] ;
         else if ( outflag == 1 )
           val = grid[i][j] / (weights[i][j] - nweights[i][j]) ;
         else if ( outflag == 2 ) val = grid[i][j] ;
         else if ( outflag == 3 ) val = grid[i][j] ; 

         /** Round and cut off neg values (if appropriate) **/

         ival = (int) ( val + CUTOFF ) ; 
         if ( cut ) 
           {
           if ( ival < 0 ) ival = 0 ;
           if ( val < 0.0 ) val = 0.0 ;
           }

         /** Track the min and max output values for each band **/ 

         if ( first == 0 ) 
           {
           min2[nb] = val ;
           max2[nb] = val ;
           first = 1 ;
           }
         else  
           {
           if ( val > max2[nb] ) max2[nb] = val ;
           if ( val < min2[nb] ) min2[nb] = val ;
           }

         }
      
       if ( datatype == 1 ) outbuf1[j] = ival ;  /* int output value   */ 
       if ( datatype == 2 ) outbuf2[j] = ival ;  /* int output value   */ 
       if ( datatype == 3 ) outbuf3[j] = ival ;  /* int output value   */ 
       if ( datatype == 4 ) outbuf4[j] = val ;   /* float output value */ 

       }

/** Write results to output file **/

     if ( datatype == 1 ) write ( opf, outbuf1, (size_t) nsout*datasize ) ;
     if ( datatype == 2 ) write ( opf, outbuf2, (size_t) nsout*datasize ) ;
     if ( datatype == 3 ) write ( opf, outbuf3, (size_t) nsout*datasize ) ;
     if ( datatype == 4 ) write ( opf, outbuf4, (size_t) nsout*datasize ) ;

     }

   if ( cut ) if ( min2[nb] < 0.0 ) min2[nb] = 0.0 ;
   printf ("\n Input Range for band %3d:  %f to %f\n", nb, min[nb], max[nb] );
   printf (" Output Range for band %3d: %f to %f\n", nb, min2[nb], max2[nb] );
   printf (" Maximum weight = %f\n",wtmax );

   }

 return ;

 }