// dem1.cpp : Defines the entry point for the console application.
//


/*
      convert dem to x, y, z format, one dem profile at a time,
      converts 3 arc second to lat/long,  sol katz, mar. 94,
      added sampling and cutting to size, sol katz, apr 94
      changed the calculation of start position in sampling code, 
                                          sol katz, jan 95
      ver 4, got rid of last 160 bytes on header line. sol katz, mar 97
*/

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

void	getheader(FILE *);
void	processProfiles(FILE *, FILE *);
void	writeNormal(FILE *, double, double);
void	writeSubset(FILE *, double, double);
void    writeCount();

#define     YMAX     2048          /*  max length of an NCIC scan line */
#define     XMAX     2048          /*  max # of NCIC scan lines */

#define MIN(x,y)     (((x) < (y)) ? (x) : (y))
#define MAX(x,y)     (((x) > (y)) ? (x) : (y))

#define SW     0
#define NW     1
#define NE     2
#define SE     3


int             base[XMAX];         /* array of base elevations */
char            mapLabel[145];
char            dump[161];
int             DEMlevel, elevationPattern, groundSystem, groundZone;
double          projectParams[15];
int             planeUnitOfMeasure, elevUnitOfMeasure, polygonSizes;
double          groundCoords[4][2], elevBounds[2], localRotation;
int             accuracyCode;
double          spatialResolution[3];
int             profileDimension[2];
int             firstRow, lastRow;
int             wcount = 0;

double          verticalScale;      /* to stretch or shrink elevations */
double          deltaY;
char            inText[24], **junk;

double          eastMost, westMost, southMost, northMost;
int             eastMostSample,  westMostSample, 
                southMostSample, northMostSample;
int             rowCount, columnCount, r, c;
long int        cellCount = 0 ;
int		rowStr=1,rowEnd, colStr=1, colEnd,
		colInt=1, rowInt=1, outType ;

void main(argc, argv)
    char           *argv[];

{
    FILE       *inf, *outf;
    double     comp = 0.0;

    if (argc != 3) {
       (void) fprintf(stderr, "Usage: dem2xyz dem_name xyz_file\n");
       exit(1);
    }

    if ((inf = fopen(argv[1], "r")) == NULL) {
       perror(argv[1]);
       exit(1);
    }
    fprintf(stderr,
         "DEM to x,y,z ascii file, version 6 -- Sol Katz, BLM, April 1997 \n ");
    fprintf(stderr,
         "         USGS DEM does not require record delimiters \n ");

    verticalScale = 1.0;

    (void)getheader(inf);


    /* 3 second DEMs overlap next one by one row, others don't */

    if (spatialResolution[0] == 3.0) {
	 comp   = spatialResolution[0] - 1.0 ;
	 deltaY = spatialResolution[1] / 3600.  ;
    }
    else {
	 comp   = 0.;
	 deltaY = spatialResolution[1] ;
    }

    eastMost  = MAX(groundCoords[NE][0], groundCoords[SE][0]);
    westMost  = MIN(groundCoords[NW][0], groundCoords[SW][0]);
    northMost = MAX(groundCoords[NE][1], groundCoords[NW][1]);
    southMost = MIN(groundCoords[SW][1], groundCoords[SE][1]);
    eastMostSample =  ((int) (eastMost / spatialResolution[0]))
		     * (int) spatialResolution[0];
    westMostSample =  ((int) ((westMost + comp) / spatialResolution[0]))
		     * (int) spatialResolution[0] ;
    northMostSample = ((int) (northMost / spatialResolution[1]))
		     * (int) spatialResolution[1] ;
    southMostSample = ((int) ((southMost + comp) / spatialResolution[1]))
		     * (int) spatialResolution[1] ;

    if (spatialResolution[0] == 3.0) { /* then it's a 1x1 degree DEM */
       printf(" eastMostSample %10f,  westMostSample %10f\n",
		eastMostSample/3600., westMostSample/3600.);
       printf("northMostSample %10f, southMostSample %10f\n",
	       northMostSample/3600., southMostSample/3600.);
       }


    columnCount = (eastMostSample  - westMostSample)
		  / (int) spatialResolution[0] + 1;
    rowCount    = (northMostSample - southMostSample)
		  / (int) spatialResolution[1] + 1;

    if (columnCount != profileDimension[1]) {
      (void) printf("       CALCULATED column Count %d != HEADER %d, ",
	    columnCount, profileDimension[1]);
      (void) printf("will use SMALLER column Count\n");
	    columnCount =MIN( profileDimension[1], columnCount);
    }

    (void) printf("number of rows %d, number of columns %d\n",
		rowCount, columnCount);

    if (argc > 2) {
      if ((outf = fopen(argv[2], "wb")) == NULL) {
	  (void) fprintf(stderr, "can't open output file %s\n", argv[2]);
	 exit(1);
       }
    }

    colInt = 1 ;
    rowInt = 1 ;
    fprintf(stderr,
	"Enter 0 for all,  1 for samples,  2 for subset : ");
    scanf("%d",&outType);

    if ( outType == 1)  /* samples */
    {
       fprintf(stderr,"Enter column sample interval :");
       scanf("%d",&colInt);
       fprintf(stderr,"Enter row    sample interval :");
       scanf("%d",&rowInt);
       /* adjust deltaY */
       deltaY = deltaY * rowInt;
    }
    else if ( outType == 2) /* subset */
    {
       fprintf(stderr,"Enter start column (from left) :");
       scanf("%d",&colStr);
       fprintf(stderr,"Enter end   column (from left) :");
       scanf("%d",&colEnd);
       fprintf(stderr,"Enter start row    (from bottom) :");
       scanf("%d",&rowStr);
       fprintf(stderr,"Enter end   row    (from bottom) :");
       scanf("%d",&rowEnd);
     }

    fprintf(stderr,"Enter a vertical scaling factor: ");
    scanf("%lf",&verticalScale);
    fprintf(stderr,"\n");

    /* adjust the end of the file */
    if ( outType == 2 && columnCount > colEnd ) columnCount = colEnd;

    /* have to read everything */
    (void) processProfiles(inf,outf);

    fprintf(stderr,"\n  %ld lines were written to the file. \n",cellCount);
/*    writeCount(); */
    exit(0);
}

/*
 * read profiles
 */

void processProfiles(inputFile,outputFile)
    FILE           *inputFile,*outputFile;
{
    int 	    c, r, mod ;
    int             count, tempInt, lastProfile = 0;

    float           noValue = 0.0;
    int             profileID[2], profileSize[2];
    double          planCoords[2], localElevation, elevExtremea[2];

    for (c = 1; c <= columnCount; c++) 
    {
       count = fscanf(inputFile, "%6d%6d%6d%6d",
		 &profileID[0],        /* 1 */
		 &profileID[1],        /* 1 */
		 &profileSize[0],      /* 2 */
		 &profileSize[1]);     /* 2 */
       if (count != 4) 
       {
	   (void) fprintf(stderr, "\nshort read of %d on column %d\n",count, c);
           printf( "id1=%6d id2=%6d size1=%6d size2=%6d\n",
		 &profileID[0],        /* 1 */
		 &profileID[1],        /* 1 */
		 &profileSize[0],      /* 2 */
		 &profileSize[1]);     /* 2 */
           fprintf(stderr,"\n  %ld lines were written to the file. \n",cellCount);
/*	   writeCount(); */
	   exit (1);
       }
       fscanf(inputFile,"%24c", inText);
       inText[20] = 'E';
       planCoords[0]   = strtod(inText,junk);     /* 3 */
       fscanf(inputFile,"%24c", inText);
       inText[20] = 'E';
       planCoords[1]   = strtod(inText,junk);     /* 3 */
       fscanf(inputFile,"%24c", inText);
       inText[20] = 'E';
       localElevation  = strtod(inText,junk);     /* 4 */
       fscanf(inputFile,"%24c", inText);
       inText[20] = 'E';
       elevExtremea[0] = strtod(inText,junk);     /* 5 */
       fscanf(inputFile,"%24c", inText);
       inText[20] = 'E';
       elevExtremea[1] = strtod(inText,junk);     /* 5 */

/*     next 2 lines are a kludge to force the end of processing */
       if ( profileID[1]-1 != lastProfile )
       {
         fprintf(stderr,"\n  %ld lines were written to the file. \n",cellCount);
         /* writeCount(); */
	 exit(0);
       }
       lastProfile =  profileID[1];

       (void) printf("Column %d has %d rows              \r",
		       profileID[1], profileSize[0]);

       firstRow = ((int) (planCoords[1] - southMostSample))
		 / (int) spatialResolution[1];
       lastRow = firstRow + profileSize[0] - 1;

       for ( r = 0 ; r < firstRow  ; r++ ) base[r] = 0;


/*     read in all the data for this column */
       for (r = firstRow; r <= lastRow; r++ )
       {
	   count = fscanf(inputFile, "%6d",	&tempInt);
	   base[r] = tempInt;
       }

/*     if cutting out a section, adjust the rows */
       if ( outType == 2 ) /* subset */ 
       {
/*    
         if ( c >= colStr ) 
           printf (" 10  outtype c=%d, colInt=%d colStr=%d columnCount=%d\n",
          c,colInt,colStr,columnCount);
*/
/*         if ( firstRow < rowStr )*/
         {
	    firstRow = rowStr ;
	    /* adjust the ycoord */
            /* suggested change */
            planCoords[1] = planCoords[1] + ( firstRow -1 ) * (deltaY);

            /* ssk dec 3, 95:  i haven't checked if we need (deltaY / 3600.);*/
            /* old version   
            planCoords[1] = planCoords[1] + ( rowInt -1 ) * deltaY;
            */
         }
         if ( lastRow  > rowEnd ) lastRow  = rowEnd ;
/*         printf (" 2 firstrow = %d, lastrow =%d, c=%d \n",firstRow,lastRow,c);*/
       }

       mod = c % colInt;
 
       if ( mod == 0 && c >= colStr )
       {
         /*   fprintf(stderr,
              "\n profile = %d , colINt = %d, mod = %d\n",c, mod); */
         if ( outType == 2 ) /* subset */ 
	   writeSubset(outputFile, planCoords[0], planCoords[1]);
         else
	   writeNormal(outputFile, planCoords[0], planCoords[1]);
       }
    }       
}

/*
 * write all or sampled data to xyz file
 *
 */

void writeNormal(outputFile, XCoord, YCoord)
    FILE           *outputFile;
    double         XCoord,   YCoord;
{
    int             r, c ;
    double          tempFloat;

/*
   wcount++;
   fprintf(stderr,"\n in writeNormal %d %d %d %d \n",
           firstRow,lastRow,rowInt,wcount);
*/
   if (spatialResolution[1] == 3.0)
   {
	 /* 3 degrees of arc */
	 XCoord /=  3600. ;
	 YCoord /=  3600. ;
   }

   for (r = firstRow; r <= lastRow; r += rowInt) {


     /* now, scale the raw value */

      tempFloat = (float) base[r]  * verticalScale;

      /*   and finally write it out to the file   */

      fprintf(outputFile,"%f  %f  %d\n",
	      XCoord, YCoord, (int)tempFloat);

      cellCount++ ;
      /*    move up the delta y distance */
      YCoord += deltaY;
   }
}

/*
 * write subset data to xyz file
 *
 */

void writeSubset(outputFile, XCoord, YCoord)
    FILE           *outputFile;
    double         XCoord,   YCoord;
{
    int             r, c ;
    double          tempFloat;


/*
   wcount++;
   fprintf(stderr,"\n in writeSubset %d %d %d %d \n",
               firstRow,lastRow,rowInt,wcount);
*/
   if (spatialResolution[1] == 3.0)
   {
	 /* 3 degrees of arc */
	 XCoord /=  3600. ;
	 YCoord /=  3600. ;
   }

   for (r = firstRow; r <= lastRow; r += rowInt) {


      /* now, scale the raw value */

      tempFloat = (float) base[r]  * verticalScale;

      /*   and finally write it out to the file   */

/*      if ( r >= rowStr - 1) */
        fprintf(outputFile,"%f  %f  %d\n",
	      XCoord, YCoord, (int)tempFloat);
  
      cellCount++ ;
      /*    move up the delta y distance */
      YCoord += deltaY;
   }
}


void getheader(inputFile)
    FILE           *inputFile;
{
    int             count,k,l ;
    count = fscanf(inputFile, "%144c%6d%6d%6d%6d",
	     mapLabel,              /* 1 */
	     &DEMlevel,             /* 2 */
	     &elevationPattern,     /* 3 */
             &groundSystem,         /* 4 */
             &groundZone );         /* 5 */
/*    (void) printf(" 5 there were %d items read from the header\n", count);*/


    for ( k=0; k<15 ; k++) {
      fscanf(inputFile,"%24c", inText);    /* 6 */
      inText[20] = 'E';
      projectParams[k]  = strtod(inText,junk);
    }

    count = fscanf(inputFile,"%6d%6d%6d",
             &planeUnitOfMeasure,   /* 7 */
             &elevUnitOfMeasure,    /* 8 */
             &polygonSizes);        /* 9 */
/*    (void) printf("3 there were %d items read from the header\n", count); */

    for ( k=0; k < 4 ; k++)
           for ( l=0; l < 2 ; l++) {
             fscanf(inputFile,"%24c", inText);    /* 6 */
             inText[20] = 'E';
             groundCoords[k][l]  = strtod(inText,junk);
	   }

    fscanf(inputFile,"%24c", inText);
    inText[20] = 'E';
    elevBounds[0] = strtod(inText,junk);       /* 10 */
    fscanf(inputFile,"%24c", inText);
    inText[20] = 'E';
    elevBounds[1] = strtod(inText,junk);      /* 11 */
    fscanf(inputFile,"%24c", inText);
    inText[20] = 'E';
    localRotation  = strtod(inText,junk);      /* 12 */

/*    (void) printf("11 there were %d items read from the header\n", count); */
    count = fscanf(inputFile,"%1d%12le%12le%12le%6d%6d",
             &accuracyCode,             /* 13 */
	     &spatialResolution[0],     /* 14 */
             &spatialResolution[1],     /* 14 */
             &spatialResolution[2],     /* 14 */
             &profileDimension[0],      /* 15 */
             &profileDimension[1]);     /* 15 */
/*    (void) printf("6 there were %d items read from the header\n", count); */

           printf("Map label: %s\n",mapLabel);
    (void) printf("[ 2] DEMlevel:  %d\n",DEMlevel );
    (void) printf("[ 3] elevation Pattern: %d\n",elevationPattern);
    (void) printf("[ 4] planimetric reference system: %d\n", groundSystem);
    (void) printf("[ 5] ground Zone: %d\n", groundZone);
    (void) printf("[ 6] proj: %15.7lf %15.7lf %15.7lf\n",
           projectParams[0], projectParams[1], projectParams[2]);
    (void) printf("[ 6] proj: %15.7lf %15.7lf %15.7lf\n",
           projectParams[3], projectParams[4], projectParams[5]);
    (void) printf("[ 6] proj: %15.7lf %15.7lf %15.7lf\n",
           projectParams[6], projectParams[7], projectParams[8]);
    (void) printf("[ 6] proj: %15.7lf %15.7lf %15.7lf\n",
           projectParams[9], projectParams[10], projectParams[11]);
    (void) printf("[ 6] proj: %15.7lf %15.7lf %15.7lf\n",
           projectParams[12], projectParams[13], projectParams[14]);
    (void) printf("[ 7] plane Unit Of Measure: %d\n",planeUnitOfMeasure );
    (void) printf("[ 8] elevation measurement units code: %d\n",
                        elevUnitOfMeasure);
    (void) printf("[ 9] polygon Sides: %d\n",polygonSizes );
    (void) printf("[10] ground coordinates\n");
    (void) printf("     NW (%14.5lf, %14.5lf)  NE (%14.5lf, %14.5lf)\n",
        groundCoords[1][0], groundCoords[1][1],
        groundCoords[2][0], groundCoords[2][1]);
    (void) printf("     SW (%14.5lf, %14.5lf)  SE (%14.5lf, %14.5lf)\n",
        groundCoords[0][0], groundCoords[0][1],
        groundCoords[3][0], groundCoords[3][1]);
        printf("Min: %15.5lf, Max: %15.5lf, Accuracy code: %d\n",
                        elevBounds[0],elevBounds[1],accuracyCode);
        printf("Spatial Resolution: %15.5lf, %15.5lf, %15.5lf \n",
           spatialResolution[0],spatialResolution[1],spatialResolution[2]);

    (void) printf("[15] map size is %d x %d\n",
                   profileDimension[0], profileDimension[1]);
/* get rid of the last 160 bytes */
    fscanf(inputFile, "%160c",dump);
/*    printf("last 160 char\n>%s<\n", dump);*/
}