grid_io.c

Go to the documentation of this file.

static char *SCCSID = "@(#)grid_io.c    1.10 05/22/00 SFWMD Research and Planning Departments";

/*

Copyright, 1995, South Florida Water Management District

DISCLAIMER:

ANY INFORMATION, INCLUDING BUT NOT LIMITED TO SOFTWARE AND DATA,
RECEIVED FROM THE SOUTH FLORIDA WATER MANAGEMENT DISTRICT ("DISTRICT")
IN FULFILLMENT OF A PUBLIC RECORDS REQUEST IS PROVIDED "AS IS" WITHOUT
WARRANTY OF ANY KIND, AND THE DISTRICT EXPRESSLY DISCLAIMS ALL EXPRESS
AND IMPLIED WARRANTIES, INCLUDING BUT NOT LIMITED TO THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
THE DISTRICT DOES NOT WARRANT, GUARANTEE, OR MAKE ANY REPRESENTATIONS
REGARDING THE USE, OR THE RESULTS OF THE USE, OF THE INFORMATION
PROVIDED TO YOU BY THE DISTRICT IN TERMS OF CORRECTNESS, ACCURACY,
RELIABILITY, TIMELINESS OR OTHERWISE.  THE ENTIRE RISK AS TO THE
RESULTS AND PERFORMANCE OF ANY INFORMATION OBTAINED FROM THE DISTRICT
IS ENTIRELY ASSUMED BY THE RECIPIENT.

*/

/* -------------------------------------------------------------
   grid_io
   This module containes routines to read and write the grid
   for the SFWMM.  The following routines are included:

   write_grid_header - writes the grid definition header 
       in binary format to a file
   read_grid_header - reads the grid definition header in
       binary format from a file
   grid_write  - writes the grid values in binary format
   grid_read   - reads the grid values in binary format
   grid_count_snapshots  -  returns the num of snapshots

   This module should contain all necessary routines for 
   handling binary grid data from the SFWMM.  The functionality
   should include at least the following:

      (1) definition of a general binary format containing at 
      least (a) specification of the geographical grid (b)
      run identification (c) snapshot identification (date) 
      (d) snapshot of data for the entire grid

      (2) reading and writing of a snapshot from and to files

      (3) 
   ------------------------------------------------------------- */
#define GCC_COMPILER 1
#include <stdio.h>
#include <string.h>
#include <ctype.h>
/* ELMchange #include <varargs.h> no longer supported gcc, need stdarg.h instead */
#include <stdarg.h>
#include <stdlib.h>

#include "grid_io.h"

#define TOO_MANY_NODES 10000
#ifndef TRUE
#define TRUE 1
#define FALSE 0
#endif

/* ELMchange Endian Swapping Routines */

#ifdef i386
#define Flip_int16(type) ((((type) >>8) & 0x00ff) | (((type) <<8) & 0xff00))
#define Flip_int16_ip(type) (*(type)=(((*(type) >>8) & 0x00ff) | ((*(type) <<8) & 0xff00)))
#define Flip_int32(type) ((((type) >>24) & 0x000000ff) | (((type) >> 8) & 0x0000ff00) | (((type) << 8) & 0x00ff0000) | (((type) << 24) & 0xff000000) )
#define Flip_int32_ip(type) (*(type)=(((*(type) >>24) & 0x000000ff) | ((*(type) >> 8) & 0x0000ff00) | ((*(type) << 8) & 0x00ff0000) | ((*(type) << 24) & 0xff000000) ))
#define Flip_float_ip(type) (Flip_int32_ip((long *)(type)))
#else
#define Flip_int16(type) (type)
#define Flip_int16_ip(type) (*(type))
#define Flip_int32(type) (type)
#define Flip_int32_ip(type) (*(type))
#define Flip_float_ip(type) (*(type))
#endif

/* internal routines */
void grid_io_fatal(char *name, char *message);
void F77_to_C_string(char *dest, char *src, int length);
void C_to_F77_string(char *dest, char *src, int length);
int openfilef77_(int *unit, char *filename, char *access, int filename_len, int access_len);
int closefilef77_(int *unit);
FILE *getfilep_(int *unit);
/* data for fortran interface routines */
static GRID fortran_grid = {0};
GR_FORTRAN_FILE *gr_fortran_file_list = 0;

/* ------------------------------------------------------------- 
   grid_write_header 
   writes the grid definition header to a file
 ------------------------------------------------------------- */

int grid_write_header(FILE *file, GRID *grid)
{
  int array_size, errs;
  int *tmp_xstart, *tmp_xend, *tmp_cum_node_count;
  errs = 0;
  tmp_xstart = grid->config.xstart;
  tmp_xend = grid->config.xend;
  tmp_cum_node_count = grid->config.cum_node_count;

  {
    GR_HEADER tmpHeader;
    strncpy(tmpHeader.title, grid->header.title, GRID_TITLE_LENGTH);
    tmpHeader.number_of_rows = Flip_int32(grid->header.number_of_rows);
    tmpHeader.number_of_nodes = Flip_int32(grid->header.number_of_nodes);
    tmpHeader.size.x =  grid->header.size.x;
    Flip_float_ip(&(tmpHeader.size.x));
    tmpHeader.size.y =  grid->header.size.y;
    Flip_float_ip(&(tmpHeader.size.y));
    if (fwrite(&tmpHeader, sizeof(GR_HEADER), 1, file) == 0) 
      ++errs;
  }

  if (grid->header.number_of_rows < MAX_GRID_ROWS)
    array_size = MAX_GRID_ROWS;
  else
    array_size = grid->header.number_of_rows;

#ifdef i386
  {
    int i;
    tmp_xstart =  (int *) malloc(array_size * sizeof(int));
    memcpy(tmp_xstart, grid->config.xstart, array_size * sizeof(int));

    tmp_xend =  (int *) malloc(array_size * sizeof(int));
    memcpy(tmp_xend, grid->config.xend, array_size * sizeof(int));

    tmp_cum_node_count =  (int *) malloc(array_size * sizeof(int));
    memcpy(tmp_cum_node_count, grid->config.cum_node_count, array_size * sizeof(int));

    for (i=0; i<array_size; ++i) {
      Flip_int32_ip(tmp_xstart+i);
      Flip_int32_ip(tmp_xend+i);
      Flip_int32_ip(tmp_cum_node_count+i);
    }
  }
#endif
  if (fwrite(tmp_xstart, sizeof(int), array_size, file) == 0)
    ++errs;
  if (fwrite(tmp_xend, sizeof(int), array_size, file) == 0)
    ++errs;
  if (fwrite(tmp_cum_node_count, sizeof(int), array_size, file) == 0)
    ++errs;

  if (errs)
    grid_io_fatal("write_grid_header", "Unable to write file header\n");

#ifdef i386
  free(tmp_xstart); free(tmp_xend); free(tmp_cum_node_count);
#endif
  
  return (0);
}


/* -------------------------------------------------------------
   wgridhd_
   Fortran interface to the grid_write_header
   ------------------------------------------------------------- */

void gridwhd_(int *fd, int *errs)
{
  FILE *file;
  file = getfilep_(fd);   

  *errs = grid_write_header(file, &fortran_grid);
}

/* -------------------------------------------------------------
   grid_read_header
   reads grid definition information from a file
   returns 0 on success, -1 on eof, or some positive number on
   error
   ------------------------------------------------------------- */

int grid_read_header(FILE *file, GRID *grid)
{
  int errs = 0, num_read, array_size;

  if ((num_read = fread(&grid->header, sizeof(GR_HEADER), 1, file)) == 0) {
    if (feof(file))
      errs = -1;
    else 
      grid_io_fatal("read_grid_header", "Unable to read header from file\n");
  }
  Flip_int32_ip(&(grid->header.number_of_rows));
  Flip_int32_ip(&(grid->header.number_of_nodes));
  Flip_float_ip(&(grid->header.size.x));
  Flip_float_ip(&(grid->header.size.y));

  if (grid->header.number_of_rows < MAX_GRID_ROWS)
    array_size = MAX_GRID_ROWS;
  else
    array_size = grid->header.number_of_rows;

  /* allocate space for configuration arrays */
  grid->config.xstart = (int *) malloc(array_size * sizeof(int));

  grid->config.xend = (int *) malloc(array_size * sizeof(int));

  grid->config.cum_node_count = (int *) malloc(array_size * sizeof(int));

  if ((num_read = fread(grid->config.xstart, sizeof(int), array_size, file)) == 0) {
    if (feof(file))
      errs = -1;
    else 
      grid_io_fatal("read_grid_header", "Unable to read header from file\n");
  }

  if ((num_read = fread(grid->config.xend, sizeof(int), array_size, file)) == 0) {
    if (feof(file))
      errs = -1;
    else
      grid_io_fatal("read_grid_header", "Unable to read header from file\n");
  }

  if ((num_read = fread(grid->config.cum_node_count, sizeof(int), array_size, file)) == 0) {
    if (feof(file))
      errs = -1;
    else
      grid_io_fatal("read_grid_header", "Unable to read header from file\n");
  }
  
#ifdef i386
  {
    int i;
    for (i=0; i<array_size; ++i) {
      Flip_int32_ip(grid->config.xstart+i);
      Flip_int32_ip(grid->config.xend+i);
      Flip_int32_ip(grid->config.cum_node_count+i);
    }
  }
#endif
  return(errs);
}


void gridrhd_(int *fd, int *errs)
{
  FILE *file;
  file = getfilep_(fd);   
  *errs = grid_read_header(file, &fortran_grid);
}

/* -------------------------------------------------------------
   grid_write
   writes a ``snapshot'' of areal data (defined by grid) to a 
   binary file 
   ------------------------------------------------------------- */
int grid_write(FILE *file, GRID *grid, char *tag, float *values)
{
  char *char_ptr;
  int written, errs = 0;
  float *tmp_values = values;

  if ((written = fwrite(tag, sizeof(char), GRID_TAG_LENGTH, file)) == 0)
    grid_io_fatal("grid_write", "Unable to write grid tag\n");
  if (written < GRID_TAG_LENGTH) errs++;
#ifdef i386
  {
    int i;
    tmp_values = (float *) malloc(grid->header.number_of_nodes * sizeof(float));
    memcpy(tmp_values, values, grid->header.number_of_nodes * sizeof(float));
    for (i=0; i<grid->header.number_of_nodes; ++i) {
      Flip_float_ip(tmp_values+i);
    }
  }
#endif
    if ((written = fwrite(tmp_values, sizeof(float), grid->header.number_of_nodes, file)) == 0)
      grid_io_fatal("grid_write", "Unable to write grid data\n");

  if (written < grid->header.number_of_nodes) errs++;
#ifdef i386
  free (tmp_values);
#endif
  return(errs);

}

/* -------------------------------------------------------------
   gwrite_
   Fortran interface to the grid_write routine
   ------------------------------------------------------------- */
void gwrite_(int *fd, char *tag, float *values, int *errs, int tag_len)
{
  char buf[GRID_TAG_LENGTH];
  FILE *file;
  
  file = getfilep_(fd);   
  F77_to_C_string(buf, tag, ((tag_len < GRID_TAG_LENGTH) ? tag_len : GRID_TAG_LENGTH));
  *errs = grid_write(file, &fortran_grid, buf, values);
}

/* -------------------------------------------------------------
   grid_read
   reads a ``snapshot'' of areal data
   returns 0 on success, -1 on eof, and some positive number on
   partial read
   ------------------------------------------------------------- */
int grid_read(FILE *file, GRID *grid, char *tag, float *values)
{
  int num_read, errs = 0;

  /* read the tag */
  if ((num_read = fread(tag, sizeof(char), GRID_TAG_LENGTH, file)) == 0) {
    if (feof(file))
      return (-1);
    else
      grid_io_fatal("grid_read", "Uable to read grid tag\n");
  }
  if (num_read != GRID_TAG_LENGTH) errs++;

  /* read the values */
  if ((num_read = fread(values, sizeof(float), grid->header.number_of_nodes, file)) == 0) 
    grid_io_fatal("grid_read", "Unable to read data\n");
#ifdef i386
  {
    int i;
    for (i=0; i<grid->header.number_of_nodes; ++i) {
      Flip_float_ip(values+i);
    }
  }
#endif
  if (num_read != grid->header.number_of_nodes) errs++;
  return (errs);
}

/* -------------------------------------------------------------
   gread_
   Fortran interface to the grid_read routine
   ------------------------------------------------------------- */
void gread_(int *fd, char *tag, float *values, int *errs, int tag_len)
{
  char buf[GRID_TAG_LENGTH];
  FILE *file;
  
  file = getfilep_(fd);   
  *errs = grid_read(file, &fortran_grid, buf, values);
  C_to_F77_string(tag, buf, ((tag_len < GRID_TAG_LENGTH) ? tag_len : GRID_TAG_LENGTH));
}

/* -------------------------------------------------------------
   grid_skip
   This routine will move the file pointer the specified number
   of records.
   ------------------------------------------------------------- */
int grid_skip(FILE *file, GRID *grid, int count)
{
  int header_size, array_size;
  int errs = 0;
  long int rec_len =  GRID_TAG_LENGTH*sizeof(char) +
    grid->header.number_of_nodes*sizeof(float);
  long int end, current;

  if (count < 0) {
    if (grid->header.number_of_rows < MAX_GRID_ROWS)
      array_size = MAX_GRID_ROWS;
    else
      array_size = grid->header.number_of_rows;
    header_size = sizeof(GRID) + 3 * sizeof(int) * array_size;
    if (ftell(file) <  header_size + abs(count)*rec_len)
      errs = grid_top(file, grid);
    else
      fseek(file, count*rec_len, 1);
  } else if (count > 0) {
    current = ftell(file);
    fseek(file, 0, 2);
    end = ftell(file);
    if ((end - current)/rec_len < count)
      errs = grid_bottom(file, grid);
    else 
      errs = fseek(file, (long)(current + count*rec_len), 0);
  }
  
  return(errs);
}

/* -------------------------------------------------------------
   gridskp_
   Fortran interface to the grid_skip routine.  Note
   that it changes the contents of the fortran_grid.  To get
   the header info a call to getgrid_ is needed (cjn 1/94)
   ------------------------------------------------------------- */
void gridskp_(int *fd, int *cnt, int *errs)
{
  FILE *file; 
  file = getfilep_(fd);   
  *errs = grid_skip(file, &fortran_grid, *cnt);
}

/* -------------------------------------------------------------
   grid_top
   This routine places the file pointer before the first data
   record in the file
   ------------------------------------------------------------- */
int grid_top(FILE *file, GRID *grid)
{
  fseek(file, 0L, 0);
  return(grid_read_header(file, grid));
}

/* -------------------------------------------------------------
   grid_bottom
   This routine moves the file pointer to just before the final
   data record in the file.
   ------------------------------------------------------------- */
int grid_bottom(FILE *file,GRID *grid)
{
  fseek(file, 0L, 2);
  return(grid_skip(file, grid, -1));
}

/* -------------------------------------------------------------
   grid_count_snapshots
   This routine returns the num of snapshots in the given gridio binary
   file. It assumes that the gridio binary file has already been opened
   and the grid_header already read. Before returning, it sets the
   file pointer back to the 1st record.
   ------------------------------------------------------------- */
int grid_count_snapshots(FILE *file, GRID *grid) {
  int i = 0;
  char tag[GRID_TAG_LENGTH];
  float *data = (float *)malloc(grid->header.number_of_nodes*sizeof(float));

  grid_top(file, grid);
  while (grid_read(file, grid, tag, data) == 0)
    i++;
  free(data);
  grid_top(file, grid);
  return(i);
}

char *strsed(register char *string, register char *pattern, int *range);

/* ELMchange (?) this grid_tag_search function was entirely commented out (?done so in original code?) */
/* -------------------------------------------------------------
   grid_tag_search
   This routine searches the grid_tags from the current position
   in the file for the regular expression passed.  The routine
   returns 1 if a matching tag was found, 0 if not.  If the
   the search was unsuccessfull the file is set to the bottom
   grid data.
   ------------------------------------------------------------- */
/*
int grid_tag_search(FILE *file, GRID *grid, char *string)
{
  char *char_ptr, tag[GRID_TAG_LENGTH], buffer[GRID_TAG_LENGTH];
  int range[2];
  float *data = (float *)malloc(grid->header.number_of_nodes*sizeof(float));
  int found = FALSE;

  sprintf(buffer, "/%s/", string);

  while (!found) {
    if (grid_read(file, grid, tag, data) != 0) {
      grid_bottom(file, grid);
      break;
    } else {
      if (strsed(tag, buffer, range) == 0)
        grid_io_fatal("grid_tag_search", "Error using strsed routine");
      else if (range[0] != -1 && range[1] != -1) {
         found = TRUE;
         grid_skip(file, grid, -1);
       }
    }
  }

  char_ptr = (char *) data;
  free(char_ptr);
  return (found);
} */


/* -------------------------------------------------------------
   grid_node
   This routine returns the array index of the node cooresponding
   to the row and column number passed.
   ------------------------------------------------------------- */
int grid_node(GRID *grid, int row, int column)
{
  if (row > 0 && row <= grid->header.number_of_rows) {
    if (column >= grid->config.xstart[row - 1] && column <= grid->config.xend[row - 1]) 
      return ((grid->config.cum_node_count[row - 1] - 1) + column - grid->config.xstart[row - 1] + 1);
  }
  return -1;
}

/* -------------------------------------------------------------
   setgrid_
   Fortran call to set internal grid definition record 
   ------------------------------------------------------------- */
void setgrid_(char *title, int *nrows, int *nnodes, float *xsize, float *ysize, int *xstart, int *xend, int *cum_count, int title_len)
{
  int i, array_size;

  F77_to_C_string(fortran_grid.header.title, title, 
                  ((title_len < GRID_TITLE_LENGTH) ? title_len : GRID_TITLE_LENGTH));

  fortran_grid.header.number_of_rows = *nrows;
  fortran_grid.header.number_of_nodes = *nnodes;
  fortran_grid.header.size.x = *xsize;
  fortran_grid.header.size.y = *ysize;

  /* allocate space for configuration arrays */
  array_size = (*nrows < MAX_GRID_ROWS) ? MAX_GRID_ROWS : *nrows;
  fortran_grid.config.xstart = (int *) malloc(array_size * sizeof(int));
  fortran_grid.config.xend = (int *) malloc(array_size * sizeof(int));
  fortran_grid.config.cum_node_count = (int *) malloc(array_size * sizeof(int));
  
  for (i = 0; i < fortran_grid.header.number_of_rows; i++) {
    fortran_grid.config.xstart[i] = xstart[i];
    fortran_grid.config.xend[i] = xend[i];
    fortran_grid.config.cum_node_count[i] = cum_count[i];
  }
}

/* -------------------------------------------------------------
   getgrid_
   Fortran call to get internal grid definition record 
   by Karen Lythgoe and Cal Neidrauer, since Bill Perkins
   documented a non-existent routine.  Documentation is like sex,
   too much is much better than not enough!
   ------------------------------------------------------------- */
void getgrid_(char *title, int *nrows, int *nnodes, float *xsize, float *ysize,
         int *xstart, int *xend, int *cum_count, int title_len)
{
  int i;

  C_to_F77_string(title,fortran_grid.header.title,  
                  ((title_len < GRID_TITLE_LENGTH) ? title_len : GRID_TITLE_LENGTH));

  *nrows = fortran_grid.header.number_of_rows;
  *nnodes = fortran_grid.header.number_of_nodes;
  *xsize = fortran_grid.header.size.x;
  *ysize = fortran_grid.header.size.y;
  for (i = 0; i < fortran_grid.header.number_of_rows; i++) {
    xstart[i] = fortran_grid.config.xstart[i];
    xend[i] = fortran_grid.config.xend[i];
    cum_count[i] =fortran_grid.config.cum_node_count[i];
  }
}

/* -------------------------------------------------------------
   grid_io_fatal
   prints fatal error messages
   ------------------------------------------------------------- */
void grid_io_fatal(char *name, char *message)
{ 
  (void) fprintf(stderr, "ERROR in %s: ", name);
  (void) fprintf(stderr, message);
  (void) abort();
}

/* -------------------------------------------------------------
   F77_to_C_string
   converts fortran style strings to C style strings
   ------------------------------------------------------------- */
void F77_to_C_string(char *dest, char *src, int length)
{
  int i;
  strncpy(dest, src, length);
  for (i = length - 1; isspace(dest[i]); i--)
    dest[i] = '\0';
}


/* -------------------------------------------------------------
   C_to_F77_string
   converts C style strings to Fortran style
   ------------------------------------------------------------- */
void C_to_F77_string(char *dest, char *src, int length)
{
  int i;

  strncpy(dest,src,length);
  for (i = strlen(dest); i < length; i++)
    dest[i] = ' ';
}


/* -------------------------------------------------------------
   grid_free
   free up memory allocated for grid header
   ------------------------------------------------------------- */
int grid_free(GRID *grid)
{
  int errs = 0;

  free(grid->config.xstart);
  free(grid->config.xend);
  free(grid->config.cum_node_count);

  return(errs);
}

int openfilef77_(int *unit, char *filename, char *access, int filename_len, int access_len)
{
  GR_FORTRAN_FILE *iFile=gr_fortran_file_list;
  GR_FORTRAN_FILE *prevFile=gr_fortran_file_list;
  char *filename_buf = (char *) malloc(sizeof(char)*(filename_len+1));
  char *access_buf = (char *) malloc(sizeof(char)*(access_len+1));

  F77_to_C_string(filename_buf, filename, filename_len);
  F77_to_C_string(access_buf, access, access_len);

  while (iFile != 0) {
    if(iFile->unit_number == *unit) {
       fprintf(stderr,"Error: unit %d already open, can not open file %s as unit %d\n", *unit, filename, *unit);
      return 0;

  }
    prevFile=iFile;
    iFile = iFile->next;
  }

  iFile = (GR_FORTRAN_FILE *) malloc(sizeof(GR_FORTRAN_FILE));
  if(prevFile == 0 )
    gr_fortran_file_list = iFile;
  else
    prevFile->next = iFile;
  iFile->fptr = fopen(filename_buf, access_buf);
  if (iFile->fptr == 0) {
     fprintf(stderr,"Error: can not open file \"%s\" as unit %d\n",
             filename_buf, *unit);
  }
  iFile->unit_number = *unit;
  iFile->next = 0;
  return 1;
}

int closefilef77_(int *unit)
{
  GR_FORTRAN_FILE *iFile=gr_fortran_file_list;
  GR_FORTRAN_FILE *prevFile=gr_fortran_file_list;
  while (iFile != 0) {
    if(iFile->unit_number == *unit) {
      if(iFile == gr_fortran_file_list)
        gr_fortran_file_list = iFile->next;
      else
        prevFile->next = iFile->next;
      free(iFile);
      return 1;
    }
  }
  return 0;
}

FILE* getfilep_(int *unit)
{
  GR_FORTRAN_FILE *iFile;
  for(iFile = gr_fortran_file_list; iFile!=0; iFile =iFile->next) {
    if(iFile->unit_number == *unit) return iFile->fptr;
  }
  return 0;
}