/* Global constants used */
#define MAXLN 4096
#define RPSHRDTYPE  1
#define RPINTDTYPE  2
#define RPFLTDTYPE  3
#include <stdio.h>
#include <stdlib.h>

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);
/*	However sometimes freeing the memory causes heap problems so I comment it out.
    free(elev2[0]); 
	free(elev2);       */
}