/*  File: demex.c 
	Program to extract subset from DEM
	Begining date: 98/1/7


*/
#include "gridio.h"

void main()
{
  char file1[MAXLN],file2[MAXLN];
  int i1,i2,j1,j2,datatype,nx,ny,err,filetype,i,j;
  FILE *fp;
  double bndbox[4],csize;
  float dx,dy,mval;


  fp = fopen("demex.in","r");
  fscanf(fp,"%d",&datatype); eol(fp);
  fscanf(fp,"%s",file1); eol(fp);
  fscanf(fp,"%s",file2); eol(fp);
  fscanf(fp,"%d%d%d%d",&i1,&i2,&j1,&j2);
  i1=i1-1;
  j1=j1-1;   /*  C 0 indexing   */
 
  if (datatype == RPSHRDTYPE)
  {
	short **elev, **elev2;
	elev2 = (short **) matalloc((j2-j1), (i2-i1), datatype);
	err=gridread(file1,(void ***)&elev,datatype,&nx,&ny,&dx,&dy,
             bndbox,&csize,&mval,&filetype);
	for(i=i1; i< i2; i++)for(j=j1; j<j2; j++)elev2[j-j1][i-i1]=elev[j][i];
	/*  Need to recompute limits here   */   
	bndbox[0]=bndbox[0]+csize*j1;        /*  xll center  */
	bndbox[1]=bndbox[1]+csize*(ny-i2);   /*  yll center  */
	bndbox[2]=bndbox[0]+csize*(j2-j1-1);  /*  other corners  */
	bndbox[3]=bndbox[1]+csize*(i2-i1-1);
	err = gridwrite(file2,(void **)elev2,datatype,j2-j1,i2-i1,dx,dy,
             bndbox,csize,mval,filetype); 
	free(elev2[0]); free(elev2);
	free(elev[0]); free(elev);
  }
    else if (datatype == RPINTDTYPE)
  {
	int **elev, **elev2;
	elev2 = (int **) matalloc((j2-j1), (i2-i1), datatype);
	err=gridread(file1,(void ***)&elev,datatype,&nx,&ny,&dx,&dy,
             bndbox,&csize,&mval,&filetype);
	for(i=i1; i< i2; i++)for(j=j1; j<j2; j++)elev2[j-j1][i-i1]=elev[j][i];
	/*  Need to recompute limits here   */   
	bndbox[0]=bndbox[0]+csize*j1;        /*  xll center  */
	bndbox[1]=bndbox[1]+csize*(ny-i2);   /*  yll center  */
	bndbox[2]=bndbox[0]+csize*(j2-j1-1);  /*  other corners  */
	bndbox[3]=bndbox[1]+csize*(i2-i1-1);
	err = gridwrite(file2,(void **)elev2,datatype,j2-j1,i2-i1,dx,dy,
             bndbox,csize,mval,filetype); 
	free(elev2[0]); free(elev2);
	free(elev[0]); free(elev);
  }
	  else if (datatype == RPFLTDTYPE)
  {
	float **elev, **elev2;
	elev2 = (float **) matalloc((j2-j1), (i2-i1), datatype);
	err=gridread(file1,(void ***)&elev,datatype,&nx,&ny,&dx,&dy,
             bndbox,&csize,&mval,&filetype);
	for(i=i1; i< i2; i++)for(j=j1; j<j2; j++)elev2[j-j1][i-i1]=elev[j][i];
	/*  Need to recompute limits here   */   
	bndbox[0]=bndbox[0]+csize*j1;        /*  xll center  */
	bndbox[1]=bndbox[1]+csize*(ny-i2);   /*  yll center  */
	bndbox[2]=bndbox[0]+csize*(j2-j1-1);  /*  other corners  */
	bndbox[3]=bndbox[1]+csize*(i2-i1-1);
	err = gridwrite(file2,(void **)elev2,datatype,j2-j1,i2-i1,dx,dy,
             bndbox,csize,mval,filetype); 
	free(elev2[0]); free(elev2);
	free(elev[0]); free(elev);
  }
	  else
	  {
		  printf("Unknown datatype %d in first row of input file\n",datatype);
		  printf("Valid datatypes are:\n1:  Short\n2:  Integer\n3:  Float\n"); 
	  }
}