 /*
  * Khoros: $Id: dem.c,v 1.2 1991/12/18 09:49:23 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: dem.c,v 1.2 1991/12/18 09:49:23 dkhoros Exp $";
#endif

 /*
  * $Log: dem.c,v $
 * Revision 1.2  1991/12/18  09:49:23  dkhoros
 * HellPatch3
 *
  */ 


/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, University of New Mexico.  All rights reserved.

 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"	 /* Copyright 1990 by UNM */
#include "vinclude.h"	
#include "file_formats/dem.h"


/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>	    file name: dem.c                                  <<<<
   >>>>                                                       <<<<
   >>>>   description: reads a USGS DEM file                  <<<<
   >>>>                                                       <<<<
   >>>>      routines: readdem				      <<<<
   >>>>                                                       <<<<
   >>>> modifications:					      <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */


#define STRINGLENGTH 50
int readfixed();

/****************************************************************
*
* Routine: readdem
*
* Purpose: This function reads a USGS DEM file and returns a ponter
*          to a structure type demimage*.
*
*   Input: The name of the DEM file.
*    
*  Output: A pointer to a structure type demimage.  NULL will be
*           returned on failure.
*    
*  Author: Per Lysne  7/31/90
*    
****************************************************************/
struct demimage *
readdem(filename)
char *filename;
{
    int i, j;
    struct demimage *newimage;
    FILE *fpt, *fopen();

    /*
     * Test for a bad input file name.
     */
    if (filename == NULL) {
        fprintf (stderr, "error in readdem - bad input file name\n");
        return (NULL);
    }

    /*
     * Open the file to read the DEM image from.
     */
    fpt = fopen (filename, "r");
    if (fpt == NULL) {
        fprintf (stderr, "error in readdem - could not access %s\n", filename);
        return(NULL);
    }

    /*
     * Allocate the structure for the DEM data.
     */
    newimage = (struct demimage *) kmalloc (sizeof(struct demimage));
    if (newimage == NULL) {
        fprintf (stderr, "error in readdem - could not allocate structure\n");
        return (NULL);
    }

    /*
     * Read the A record (header record) from the DEM file.
     */
    if (rd_readA(fpt, newimage) != 1) {
        fprintf (stderr, "error in readdem - could not read A record\n");
        free (newimage);
        return (NULL);
    }

    /*
     * Allocate the array within the DEM structure with is an array of
     * pointers to the individual columns.
     */
    newimage->data = (struct demprofile **)
                      kcalloc (newimage->numcolumns, sizeof(struct dem *));
    if (newimage->data == NULL) {
        fprintf (stderr, 
                 "error in readdem - could not allocate profile array\n");
        free (newimage);
        return (NULL);
    }
 
    /*
     * This loop reads the profiles (columns) of the data itself.
     */
    for (i=0; i<newimage->numcolumns; i++) {

        /*
         * Allocate the record for a single column of data.
         */
        newimage->data[i] = (struct demprofile *)
                             kmalloc (sizeof(struct demprofile));
        if (newimage->data[i] == NULL) {
            fprintf (stderr,
                "error in readdem - could not allocate profile no. %d\n", i);
            for (j=0; j<i; j++) {
                free (newimage->data[j]);
            }
            free (newimage->data);
            free (newimage);
            return (NULL);
        }

        /*
         * Read a single column of data.
         */
        if (rd_readB(fpt, newimage->data[i]) != 1) {
            fprintf (stderr, 
                   "error in readdem - could not read B record no. %d\n", i);
            for (j=0; j<=i; j++) {
                free (newimage->data[j]);
            }
            free (newimage->data);
            free (newimage);
            return (NULL);
        }
    }

    /*
     * Return successfully.
     */
    return (newimage);
}



/****************************************************************
*
* Routine: rd_readA  (internal routine)
*
* Purpose: This function reads the A record (header) from a USGS
*          DEM file.
*
*   Input: The inputs are a file pointer to the file to read from
*          and a pointer to a struct demimage to read the data into.
*    
*  Output: The output is the data that has been read into the structure
*          that was input.  This function returns 1 upon success and
*          0 on failure.
*    
*  Author: Per Lysne  7/31/90
*    
****************************************************************/
int
rd_readA(file, image)
FILE *file;
struct demimage *image;
{
    int i, j;

    /*
     * Check for a bad input file pointer.
     */
    if (file == NULL) {
        fprintf (stderr, "error in rd_readA - bad file input\n");
        return (0);
    }

    /*
     * Check for a bad input structure pointer.
     */
    if (image == NULL) {
        fprintf (stderr, "error in rd_readA - bad image input\n");
        return (0);
    }

    /*
     * Read the name and location comment from the top of the DEM file.
     */
    if (fgets (image->name, 145, file) == NULL) {
        fprintf (stderr, "error in rd_readA - could not read name field\n");
        return (0);
    }

    /*
     * Read the level code of the data in this file.
     */
    if (readfixed(file, 6, "%d", &(image->level)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read level field\n");
        return (0);
    }

    /*
     * Read the sampling pattern.
     */
    if (readfixed(file, 6, "%d", &(image->pattern)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read pattern\n");
        return (0);
    } 
    if ((image->pattern != REGULAR) && (image->pattern != IRREGULAR)) {
        fprintf (stderr, "warning: bad level was read\n");
        return (0);
    }

    /*
     * Read the type of ground reference system used.
     */
    if (readfixed(file, 6, "%d", &(image->reference)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read reference\n");
        return (0);
    } 
    if ((image->reference != GEOGRAPHIC) && (image->reference != UTM) &&
       (image->reference != STATE_PLANE)) {
        fprintf (stderr, "warning: bad refernce was read\n");
        return (0);
    }

    /*
     * Read the zone for the data.
     */
    if (readfixed(file, 6, "%d", &(image->zone)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read zone\n");
        return (0);
    } 
    if (image->zone < 0) {
        fprintf (stderr, "warning: bad zone was read\n");
        return (0);
    }

    /*
     * Read the 15 map projection parameters.
     */
    for (i=0; i<15; i++) {
        if (readfixed(file, 24, "%lf", &(image->projection[i])) != 1) {
            fprintf (stderr,
                     "error in rd_readA - could not read projection[%d]\n", i);
            return (0);
        }
    }

    /*
     * Read the type of units used to measure ground distances.
     */
    if (readfixed(file, 6, "%d", &(image->groundunit)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read ground unit\n");
        return (0);
    }
    if ((image->groundunit != RADIANS) && (image->groundunit != FEET) &&
       (image->groundunit != METERS) && (image->groundunit != ARC_SECONDS)) {
        fprintf (stderr, "warning: bad ground unit was read\n");
        return (0);
    }

    /*
     * Read the type of units used to measure elevation distances.
     */
    if (readfixed(file, 6, "%d", &(image->elevationunit)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read elevation unit\n");
        return (0);
    }
    if ((image->elevationunit != FEET) && (image->elevationunit != METERS)) {
        fprintf (stderr, "warning: bad elevation unit was read\n");
        return (0);
    }

    /*
     * Read the number of sides on the polygon bounding the data.
     */
    if (readfixed(file, 6, "%d", &(image->sides)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read no. of sides\n");
        return (0);
    }
    if (image->sides != 4) {
        fprintf (stderr, "warning: bounding polygon should have 4 sides\n");
        return (0);
    }

    /*
     * Read the ground coordinates of the four corners of the data area.
     */
    for (i=0; i<4; i++) {
        for (j=0; j<2; j++) {
            if (readfixed(file, 24, "%lf", &(image->corners[i][j])) != 1) {
                fprintf (stderr,
                "error in rd_readA - could not read coordinate[%d][%d]\n",i,j);
                return (0);
            }
        }
    }

    /*
     * Read the minimum elevation in the data.
     */
    if (readfixed(file, 24, "%lf", &(image->minimum)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read minimum\n");
        return (0);
    }

    /*
     * Read the maximum elevation in the data.
     */
    if (readfixed(file, 24, "%lf", &(image->maximum)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read maximum\n");
        return (0);
    }

    /*
     * Read the angle between local and universal coordinates.
     */
    if (readfixed(file, 24, "%lf", &(image->angle)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read angle\n");
        return (0);
    }

    /*
     * Read the flag that tells us if accuracy data in available.
     */
    if (readfixed(file, 6, "%d", &(image->accuracy)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read accuracy\n");
        return (0);
    }
    if ((image->accuracy != KNOWN) && (image->accuracy != UNKNOWN)) {
        fprintf (stderr, "warning: bad accuracy flag was read\n");
        return (0);
    }

    /*
     * Read the x resolution.
     */
    if (readfixed(file, 12, "%lf", &(image->xres)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read xres\n");
        return (0);
    }
    if (image->xres <= 0) {
        fprintf (stderr, "warning: bad x resolution was read\n");
        return (0);
    }

    /*
     * Read the y resolution.
     */
    if (readfixed(file, 12, "%lf", &(image->yres)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read yres\n");
        return (0);
    }
    if (image->yres <= 0) {
        fprintf (stderr, "warning: bad y resolution was read\n");
        return (0);
    }

    /*
     * Read the z resolution.
     */
    if (readfixed(file, 12, "%lf", &(image->zres)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read zres\n");
        return (0);
    }
    if (image->zres <= 0) {
        fprintf (stderr, "warning: bad z resolution was read\n");
        return (0);
    }

    /*
     * Read the number of rows in a profile.
     */
    if (readfixed(file, 6, "%d", &(image->numrows)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read no. of rows\n");
        return (0);
    }
    if (image->numrows != 1) {
        fprintf (stderr, "warning: bad number of rows per profile was read\n");
        return (0);
    }

    /* 
     * Read the number of profiles in the data.
     */
    if (readfixed(file, 6, "%d", &(image->numcolumns)) != 1) {
        fprintf (stderr, "error in rd_readA - could not read no. of columns\n");
        return (0);
    }
    if (image->numcolumns < 1) {
        fprintf (stderr, "warning: bad number of columns was read\n");
        return (0);
    }

    /*
     * Return successfully.
     */
    return(1);
}



/****************************************************************
*
* Routine: rd_readB (internal routine)
*
* Purpose: This function reads a B record (profile/column) from a
*          USGS DEM file.
*
*   Input: The input is a file pointer to the file to read from,
*          and a pointer to a struct type demprofile to read the
*          data into.
*    
*  Output: The output is the data read into the structure that was
*          input.  This function returns 1 upon success and 0 on failure.
*    
*  Author: Per Lysne  7/31/90
*    
****************************************************************/
int
rd_readB(file, profile)
FILE *file;
struct demprofile *profile;
{
    int i;

    /* 
     * Check for a bad input file pointer.
     */
    if (file == NULL) {
        fprintf (stderr, "error in rd_readB - bad file input\n");
        return (0);
    }

    /*
     * Check for a bad input structure pointer.
     */
    if (profile == NULL) {
        fprintf (stderr, "error in rd_readB - bad profile input\n");
        return (0);
    }

    /*
     * Read the row number of this profile.  Should be set to 1.  Scan
     * is used because we need to eat leading spaces.
     */
    if (fscanf (file, "%1d", &(profile->row)) <= 0) {
        fprintf (stderr, "error in rd_readB - could not read row\n");
        return (0);
    }
    if (profile->row != 1) {
        fprintf (stderr, "warning: bad profile row number was read\n");
        return (0);
    }

    /*
     * Read the column number of this profile.
     */
    if (readfixed(file, 6, "%d", &(profile->column)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read column\n");
        return (0);
    }
    if (profile->column < 0) {
        fprintf (stderr, "warning: bad profile column number was read\n");
        return (0);
    }

    /* 
     * Read the number of data (rows) in this profile.
     */
    if (readfixed(file, 6, "%d", &(profile->numrows)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read numrows\n");
        return (0);
    }
    if (profile->numrows < 1) {
        fprintf (stderr, "warning: bad number of data per profile was read\n");
        return (0);
    }

    /*
     * Read the number of data bands (columns) in this profile.  Should
     * be 1.
     */
    if (readfixed(file, 6, "%d", &(profile->numcolumns)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read numcolumns\n");
        return (0);
    }
    if (profile->numcolumns != 1) {
        fprintf (stderr, "warning: bad number of rows per profile was read\n");
        return (0);
    }

    /*
     * Read the geographic x coordinate of the first data in this profile.
     */
    if (readfixed(file, 24, "%lf", &(profile->xground)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read xground\n");
        return (0);
    }

    /*
     * Read the geographic y coordinate of the first data in this profile.
     */
    if (readfixed(file, 24, "%lf", &(profile->yground)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read yground\n");
        return (0);
    }

    /*
     * Read the base elevation of the data in this profile.
     */
    if (readfixed(file, 24, "%lf", &(profile->elevation)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read elevation\n");
        return (0);
    }

    /*
     * Read the minimum elevation in this profile.
     */
    if (readfixed(file, 24, "%lf", &(profile->minimum)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read minimum\n");
        return (0);
    }

    /*
     * Read the maximum elevation in this profile.
     */
    if (readfixed(file, 24, "%lf", &(profile->maximum)) != 1) {
        fprintf (stderr, "error in rd_readB - could not read maximum\n");
        return (0);
    }

    /*
     * Allocate an array to store the data.
     */
    profile->data = (int *) kcalloc (profile->numrows, sizeof(int));
    if (profile->data == NULL) {
        fprintf (stderr, "error in rd_readB - malloc for data failed\n");
        return (0);
    }

    /*
     * Read the elevation data.  Scan is used because the data occasionally
     * has a blank field in it.
     */
    for (i=0; i<profile->numrows; i++) {
        if (fscanf(file, "%d", &(profile->data[i])) <= 0) {
            fprintf (stderr,
                    "error in rd_readB - could not read data point %d\n", i);
            free (profile->data);
            return (0);
        }
    }

    /*
     * Return successfully.
     */
    return (1);
}
