#include <stdio.h>
#include <stdlib.h>

#include "gridio.h"   /*  This is the header file for the GRIDIO functions that 
                           access ESRI grids  */
/*  This program needs to be linked with gridio.c that is part of TARDEM and avgridio.lib 
    that is part of ESRI gioapi interface (a part of Spatial Analyst that is not installed
	by default so you may need to reinstall from the distribution CD).

   Alternatively you can link against gridio_null.c that is part of TARDEM and allows the 
   grid functions to work with ASCII files only, not the ESRI binary grids directly  */ 

/* The followind definitions from gridio.h are used.

#define MAXLN 4096
#define RPSHRDTYPE  1
#define RPINTDTYPE  2
#define RPFLTDTYPE  3

void **matalloc(int nx,int ny,int datatype);

int gridwrite(char *file, void **data, int datatype, int nx, int ny, float dx, 
 float dy, double bndbox[4], double csize, float ndv, int filetype);

int gridread(char *file, void ***data, int datatype, int *nx, int *ny,
 float *dx, float *dy, double bndbox[4], double *csize, float *ndv, int *filetype);

  */


void main(int argc,char **argv)
{
  char demfile[MAXLN];
  int err, nx, ny,filetype, icount,i,j;
  float **elev, **elev2, dx, dy, xndv;
  double bndbox[4],csize;

   if(argc != 2)
    {  
       printf("Usage:\n %s filename\n",argv[0]);
       printf("\n");
       printf("A command line argument that is the name of the DEM to work with is required\n");
       exit(0);  
    }
    sprintf(demfile,"%s",argv[1]);
    err=gridread(demfile,(void ***)&elev,RPFLTDTYPE,&nx,&ny,&dx,&dy,
             bndbox,&csize,&xndv,&filetype);
/*   This reads in the DEM.   This also allocates memory for the DEM. Now you can work with it
     The DEM I first got had some missing data in a stripe down the middle 
	 that I fixed by averaging  */
     elev2 = (float **) matalloc(nx, ny, RPFLTDTYPE);
/*  The statement above allocates memory the size of the grid just read.  
    This is not necessary in this case where the elev array could be overwritten
	but is given here to show how it is done in cases where it is needed.  */
	icount=0;
	for(i = 0; i<ny; i++)
	  for(j=1; j< (nx-1); j++)
	    if(elev[j][i] <= 0.)
		{
/*  Note that array subscripting is the other way round to comply with C row order
    and Fortran column order defaults.  Just remember that the first subscript refers
	to the column.  Also subscripts go 0 .. n-1 rather than 1 .. n   */
			if(elev[j-1][i] > 0. && elev[j+1][i] > 0.)
			{
	          elev2[j][i] = (float)(elev[j-1][i]+elev[j+1][i])*.5;
	          icount=icount+1;
			}
 	        else
	          elev2[j][i]=xndv;  /* Change to a no data value   */
		}
		else
			elev2[j][i]=elev[j][i];
/*   The above omits the first and last columns  */
	for(i = 0; i<ny; i++)
	{
		elev2[0][i] = elev[0][i];
		elev2[nx-1][i] = elev[nx-1][i];
	}
	printf("%d values fixed\n",icount);
	printf("Input file name for fixed DEM output\n");
	scanf("%s",&demfile);
    err = gridwrite(demfile,(void **)elev2,RPFLTDTYPE,nx,ny,dx,dy,
             bndbox,csize,xndv,filetype); 
/*  Good programming practice is to now tidy up - release memory  */
	free(elev[0]); 
	free(elev);
    free(elev2[0]); 
	free(elev2);       
}