src/dem.c

#include <stdio.h>
#include <string.h>
#include <glib.h>
#include <math.h>
#include <stdlib.h>
#include <zlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/mman.h>


#include "dem.h"
#include "file.h"

#define DEM_BLOCK_SIZE 1024
#define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))

static gboolean get_double_and_continue ( gchar **buffer, gdouble *tmp, gboolean warn )
{
  gchar *endptr;
  *tmp = g_strtod(*buffer, &endptr);
  if ( endptr == NULL|| endptr == *buffer ) {
    if ( warn )
      g_warning("Invalid DEM");
    return FALSE;
  }
  *buffer=endptr;
  return TRUE;
}


static gboolean get_int_and_continue ( gchar **buffer, gint *tmp, gboolean warn )
{
  gchar *endptr;
  *tmp = strtol(*buffer, &endptr, 10);
  if ( endptr == NULL|| endptr == *buffer ) {
    if ( warn )
      g_warning("Invalid DEM");
    return FALSE;
  }
  *buffer=endptr;
  return TRUE;
}

static gboolean dem_parse_header ( gchar *buffer, VikDEM *dem )
{
  gdouble val;
  gint int_val;
  guint i;
  gchar *tmp = buffer;

  /* incomplete header */
  if ( strlen(buffer) != 1024 )
    return FALSE;

  /* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
  while (*tmp) {
    if ( *tmp=='D')
      *tmp='E';
    tmp++;
  }

  /* skip name */
  buffer += 149;

  /* "DEM level code, pattern code, palaimetric reference system code" -- skip */
  get_int_and_continue(&buffer, &int_val, TRUE);
  get_int_and_continue(&buffer, &int_val, TRUE);
  get_int_and_continue(&buffer, &int_val, TRUE);

  /* zone */
  get_int_and_continue(&buffer, &int_val, TRUE);
  dem->utm_zone = int_val;
  /* TODO -- southern or northern hemisphere?! */
  dem->utm_letter = 'N';

  /* skip numbers 5-19  */
  for ( i = 0; i < 15; i++ ) {
    if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
      g_warning ("Invalid DEM header");
      return FALSE;
    }
  }

  /* number 20 -- horizontal unit code (utm/ll) */
  get_double_and_continue(&buffer, &val, TRUE);
  dem->horiz_units = val;
  get_double_and_continue(&buffer, &val, TRUE);
  /* dem->orig_vert_units = val; now done below */

  /* TODO: do this for real. these are only for 1:24k and 1:250k USGS */
  if ( dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS ) {
    dem->east_scale = 10.0; /* meters */
    dem->north_scale = 10.0;
    dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
  } else {
    dem->east_scale = 3.0; /* arcseconds */
    dem->north_scale = 3.0;
    dem->orig_vert_units = VIK_DEM_VERT_METERS;
  }

  /* skip next */
  get_double_and_continue(&buffer, &val, TRUE);

  /* now we get the four corner points. record the min and max. */
  get_double_and_continue(&buffer, &val, TRUE);
  dem->min_east = dem->max_east = val;
  get_double_and_continue(&buffer, &val, TRUE);
  dem->min_north = dem->max_north = val;

  for ( i = 0; i < 3; i++ ) {
    get_double_and_continue(&buffer, &val, TRUE);
    if ( val < dem->min_east ) dem->min_east = val;
    if ( val > dem->max_east ) dem->max_east = val;
    get_double_and_continue(&buffer, &val, TRUE);
    if ( val < dem->min_north ) dem->min_north = val;
    if ( val > dem->max_north ) dem->max_north = val;
  }

  return TRUE;
}

static void dem_parse_block_as_cont ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
{
  gint tmp;
  while ( *cur_row < GET_COLUMN(dem, *cur_column)->n_points ) {
    if ( get_int_and_continue(&buffer, &tmp,FALSE) ) {
      if ( dem->orig_vert_units == VIK_DEM_VERT_DECIMETERS )
        GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) (tmp / 10);
      else
        GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) tmp;
    } else
      return;
    (*cur_row)++;
  }
  *cur_row = -1; /* expecting new column */
}

static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
{
  guint n_rows;
  gint i;
  gdouble east_west, south;
  gdouble tmp;

  /* 1 x n_rows 1 east_west south x x x DATA */

  if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
    g_warning("Incorrect DEM Class B record: expected 1");
    return;
  }

  /* don't need this */
  if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;

  /* n_rows */
  if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
    return;
  n_rows = (guint) tmp;

  if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
    g_warning("Incorrect DEM Class B record: expected 1");
    return;
  }

  if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
    return;
  if ( !get_double_and_continue(&buffer, &south, TRUE) )
    return;

  /* next three things we don't need */
  if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
  if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
  if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;


  dem->n_columns ++;
  (*cur_column) ++;

  /* empty spaces for things before that were skipped */
  (*cur_row) = (south - dem->min_north) / dem->north_scale;
  if ( south > dem->max_north || (*cur_row) < 0 )
    (*cur_row) = 0;

  n_rows += *cur_row;

  g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
  GET_COLUMN(dem,*cur_column)->east_west = east_west;
  GET_COLUMN(dem,*cur_column)->south = south;
  GET_COLUMN(dem,*cur_column)->n_points = n_rows;
  GET_COLUMN(dem,*cur_column)->points = g_malloc(sizeof(gint16)*n_rows);

  /* no information for things before that */
  for ( i = 0; i < (*cur_row); i++ )
    GET_COLUMN(dem,*cur_column)->points[i] = VIK_DEM_INVALID_ELEVATION;

  /* now just continue */
  dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );


}

static void dem_parse_block ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
{
  /* if haven't read anything or have read all items in a columns and are expecting a new column */
  if ( *cur_column == -1 || *cur_row == -1 ) {
    dem_parse_block_as_header(buffer, dem, cur_column, cur_row);
  } else {
    dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
  }
}

/* return size of unzip data or 0 if failed */
/* can be made generic to uncompress zip, gzip, bzip2 data */
static guint uncompress_data(void *uncompressed_buffer, guint uncompressed_size, void *compressed_data, guint compressed_size)
{
	z_stream stream;
	int err;

	stream.next_in = compressed_data;
	stream.avail_in = compressed_size;
	stream.next_out = uncompressed_buffer;
	stream.avail_out = uncompressed_size;
	stream.zalloc = (alloc_func)0;
	stream.zfree = (free_func)0;
	stream.opaque = (voidpf)0;

	/* negative windowBits to inflateInit2 means "no header" */
	if ((err = inflateInit2(&stream, -MAX_WBITS)) != Z_OK) {
		g_warning("%s(): inflateInit2 failed", __PRETTY_FUNCTION__);
		return 0;
	}

	err = inflate(&stream, Z_FINISH);
	if ((err != Z_OK) && (err != Z_STREAM_END) && stream.msg) {
		g_warning("%s() inflate failed err=%d \"%s\"", __PRETTY_FUNCTION__, err, stream.msg == NULL ? "unknown" : stream.msg);
		inflateEnd(&stream);
		return 0;
	}

	inflateEnd(&stream);
	return(stream.total_out);
}

static void *unzip_hgt_file(gchar *zip_file, gulong *unzip_size)
{
	void *unzip_data = NULL;
	gchar *zip_data;
	struct _lfh {
		guint32 sig;
		guint16 extract_version;
		guint16 flags;
		guint16 comp_method;
		guint16 time;
		guint16 date;
		guint32 crc_32;
		guint32 compressed_size;
		guint32 uncompressed_size;
		guint16 filename_len;
		guint16 extra_field_len;
	}  __attribute__ ((__packed__)) *local_file_header = NULL;


	local_file_header = (struct _lfh *) zip_file;
	if (local_file_header->sig != 0x04034b50) {
		g_warning("%s(): wrong format\n", __PRETTY_FUNCTION__);
		g_free(unzip_data);
		goto end;
	}

	zip_data = zip_file + sizeof(struct _lfh) + local_file_header->filename_len + local_file_header->extra_field_len;
	unzip_data = g_malloc(local_file_header->uncompressed_size);
	gulong uncompressed_size = local_file_header->uncompressed_size;

	if (!(*unzip_size = uncompress_data(unzip_data, uncompressed_size, zip_data, local_file_header->compressed_size))) {
		g_free(unzip_data);
		unzip_data = NULL;
		goto end;
	}

end:
	return(unzip_data);
}

static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zip)
{
  gint i, j;
  VikDEM *dem;
  struct stat stat;
  off_t file_size;
  gint16 *dem_mem = NULL;
  gint16 *dem_file = NULL;
  const gint num_rows_3sec = 1201;
  const gint num_rows_1sec = 3601;
  gint num_rows;
  int fd = fileno(f);
  gint arcsec;

  dem = g_malloc(sizeof(VikDEM));

  dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
  dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;

  /* TODO */
  dem->min_north = atoi(basename+1) * 3600;
  dem->min_east = atoi(basename+4) * 3600;
  if ( basename[0] == 'S' )
    dem->min_north = - dem->min_north;
  if ( basename[3] == 'W' )
    dem->min_east = - dem->min_east;

  dem->max_north = 3600 + dem->min_north;
  dem->max_east = 3600 + dem->min_east;

  dem->columns = g_ptr_array_new();
  dem->n_columns = 0;

  if (fstat(fd, &stat) == -1)
    g_error("%s(): fstat failed on %s\n", __PRETTY_FUNCTION__, basename);
  if ((dem_file = mmap(NULL, stat.st_size, PROT_READ, MAP_PRIVATE, fd, 0)) == (void *) -1)
    g_error("%s(): mmap failed on %s\n", __PRETTY_FUNCTION__, basename);

  file_size = stat.st_size;
  dem_mem = dem_file;
  if (zip) {
    void *unzip_mem = NULL;
    gulong ucsize;

    if ((unzip_mem = unzip_hgt_file((gchar *)dem_file, &ucsize)) == NULL) {
      munmap(dem_file, file_size);
      g_ptr_array_free(dem->columns, TRUE);
      g_free(dem);
      return NULL;
    }

    dem_mem = unzip_mem;
    file_size = ucsize;
  }

  if (file_size == (num_rows_3sec * num_rows_3sec * sizeof(gint16)))
    arcsec = 3;
  else if (file_size == (num_rows_1sec * num_rows_1sec * sizeof(gint16)))
    arcsec = 1;
  else {
    g_warning("%s(): file %s does not have right size\n", __PRETTY_FUNCTION__, basename);
    munmap(dem_file, file_size);
    g_ptr_array_free(dem->columns, TRUE);
    g_free(dem);
    return NULL;
  }

  num_rows = (arcsec == 3) ? num_rows_3sec : num_rows_1sec;
  dem->east_scale = dem->north_scale = arcsec;

  for ( i = 0; i < num_rows; i++ ) {
    dem->n_columns++;
    g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
    GET_COLUMN(dem,i)->east_west = dem->min_east + arcsec*i;
    GET_COLUMN(dem,i)->south = dem->min_north;
    GET_COLUMN(dem,i)->n_points = num_rows;
    GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*num_rows);
  }

  int ent = 0;
  for ( i = (num_rows - 1); i >= 0; i-- ) {
    for ( j = 0; j < num_rows; j++ ) {
      GET_COLUMN(dem,j)->points[i] = GINT16_FROM_BE(dem_mem[ent]);
      ent++;
    }

  }

  if (zip)
    g_free(dem_mem);
  munmap(dem_file, stat.st_size);
  return dem;
}

#define IS_SRTM_HGT(fn) (strlen((fn))==11 && (fn)[7]=='.' && (fn)[8]=='h' && (fn)[9]=='g' && (fn)[10]=='t' && ((fn)[0]=='N' || (fn)[0]=='S') && ((fn)[3]=='E' || (fn)[3]=='W'))

VikDEM *vik_dem_new_from_file(const gchar *file)
{
  FILE *f;
  VikDEM *rv;
  gchar buffer[DEM_BLOCK_SIZE+1];

  /* use to record state for dem_parse_block */
  gint cur_column = -1;
  gint cur_row = -1;
  const gchar *basename = a_file_basename(file);

      /* FILE IO */
  f = fopen(file, "r");
  if ( !f )
    return NULL;

  if ( (strlen(basename)==11 || ((strlen(basename) == 15) && (basename[11] == '.' && basename[12] == 'z' && basename[13] == 'i' && basename[14] == 'p'))) &&
       basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
       (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W')) {
    gboolean is_zip_file = (strlen(basename) == 15);
    rv = vik_dem_read_srtm_hgt(f, basename, is_zip_file);
    fclose(f);
    return(rv);
  }

      /* Create Structure */
  rv = g_malloc(sizeof(VikDEM));

      /* Header */
  buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
  if ( ! dem_parse_header ( buffer, rv ) ) {
    g_free ( rv );
    return NULL;
  }
  /* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */

  rv->columns = g_ptr_array_new();
  rv->n_columns = 0;

      /* Column -- Data */
  while (! feof(f) ) {
     gchar *tmp;

     /* read block */
     buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';

     /* Fix Fortran-style exponentiation */
     tmp = buffer;
     while (*tmp) {
       if ( *tmp=='D')
         *tmp='E';
       tmp++;
     }

     dem_parse_block(buffer, rv, &cur_column, &cur_row);
  }

     /* TODO - class C records (right now says 'Invalid' and dies) */

  fclose(f);

  /* 24k scale */
  if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
    rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;

  /* FIXME bug in 10m DEM's */
  if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
    rv->min_east -= 100;
    rv->min_north += 200;
  }


  return rv;
}

void vik_dem_free ( VikDEM *dem )
{
  guint i;
  for ( i = 0; i < dem->n_columns; i++)
    g_free ( GET_COLUMN(dem, i)->points );
  g_ptr_array_free ( dem->columns, TRUE );
  g_free ( dem );
}

gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
{
  if ( col < dem->n_columns )
    if ( row < GET_COLUMN(dem, col)->n_points )
      return GET_COLUMN(dem, col)->points[row];
  return VIK_DEM_INVALID_ELEVATION;
}

gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
{
  gint col, row;

  if ( east > dem->max_east || east < dem->min_east ||
      north > dem->max_north || north < dem->min_north )
    return VIK_DEM_INVALID_ELEVATION;

  col = (gint) floor((east - dem->min_east) / dem->east_scale);
  row = (gint) floor((north - dem->min_north) / dem->north_scale);

  return vik_dem_get_xy ( dem, col, row );
}

static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
    gdouble east, gdouble north, /* in seconds */
    gint16 *elevs, gint16 *dists)
{
  int i;
  int cols[4], rows[4];
  struct LatLon ll[4];
  struct LatLon pos;

  if ( east > dem->max_east || east < dem->min_east ||
      north > dem->max_north || north < dem->min_north )
    return FALSE;  /* got nothing */

  pos.lon = east/3600;
  pos.lat = north/3600;

  /* order of the data: sw, nw, ne, se */
  cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);  /* sw */
  rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
  ll[0].lon = (east - fabs(fmod(east, dem->east_scale)))/3600;
  ll[0].lat = (north - fabs(fmod(north, dem->north_scale)))/3600;

  cols[1] = cols[0];         /*nw*/
  rows[1] = rows[0] + 1;
  ll[1].lon = ll[0].lon;
  ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;

  cols[2] = cols[0] + 1;     /*ne*/
  rows[2] = rows[0] + 1;
  ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
  ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600;

  cols[3] = cols[0] + 1;     /*se*/
  rows[3] = rows[0];
  ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
  ll[3].lat = ll[0].lat;

  for (i = 0; i < 4; i++) {
    if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION)
      return FALSE;
    dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
  }

  return TRUE;  /* all OK */
}

gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
{
  int i;
  gint16 elevs[4], dists[4];

  if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
    return VIK_DEM_INVALID_ELEVATION;

  for (i = 0; i < 4; i++) {
    if (dists[i] < 1) {
      return(elevs[i]);
    }
  }

  gdouble t = (gdouble)elevs[0]/dists[0] + (gdouble)elevs[1]/dists[1] + (gdouble)elevs[2]/dists[2] + (gdouble)elevs[3]/dists[3];
  gdouble b = 1.0/dists[0] + 1.0/dists[1] + 1.0/dists[2] + 1.0/dists[3];

  return(t/b);
}

gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
{
  return vik_dem_get_simple_interpol(dem, east, north);
}

void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )
{
  *col = (gint) floor((east - dem->min_east) / dem->east_scale);
  *row = (gint) floor((north - dem->min_north) / dem->north_scale);
  if ( *col < 0 ) *col = 0;
  if ( *row < 0 ) *row = 0;
}
Commit	Line	Data
	1	#include <stdio.h>
	2	#include <string.h>
	3	#include <glib.h>
	4	#include <math.h>
	5	#include <stdlib.h>
	6	#include <zlib.h>
	7	#include <sys/types.h>
	8	#include <sys/stat.h>
	9	#include <unistd.h>
	10	#include <sys/mman.h>
	11
	12
	13	#include "dem.h"
	14	#include "file.h"
	15
	16	#define DEM_BLOCK_SIZE 1024
	17	#define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))
	18
	19	static gboolean get_double_and_continue ( gchar *buffer, gdouble tmp, gboolean warn )
	20	{
	21	gchar *endptr;
	22	tmp = g_strtod(buffer, &endptr);
	23	if ( endptr == NULL\|\| endptr == *buffer ) {
	24	if ( warn )
	25	g_warning("Invalid DEM");
	26	return FALSE;
	27	}
	28	*buffer=endptr;
	29	return TRUE;
	30	}
	31
	32
	33	static gboolean get_int_and_continue ( gchar *buffer, gint tmp, gboolean warn )
	34	{
	35	gchar *endptr;
	36	tmp = strtol(buffer, &endptr, 10);
	37	if ( endptr == NULL\|\| endptr == *buffer ) {
	38	if ( warn )
	39	g_warning("Invalid DEM");
	40	return FALSE;
	41	}
	42	*buffer=endptr;
	43	return TRUE;
	44	}
	45
	46	static gboolean dem_parse_header ( gchar buffer, VikDEM dem )
	47	{
	48	gdouble val;
	49	gint int_val;
	50	guint i;
	51	gchar *tmp = buffer;
	52
	53	/* incomplete header */
	54	if ( strlen(buffer) != 1024 )
	55	return FALSE;
	56
	57	/* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
	58	while (*tmp) {
	59	if ( *tmp=='D')
	60	*tmp='E';
	61	tmp++;
	62	}
	63
	64	/* skip name */
	65	buffer += 149;
	66
	67	/* "DEM level code, pattern code, palaimetric reference system code" -- skip */
	68	get_int_and_continue(&buffer, &int_val, TRUE);
	69	get_int_and_continue(&buffer, &int_val, TRUE);
	70	get_int_and_continue(&buffer, &int_val, TRUE);
	71
	72	/* zone */
	73	get_int_and_continue(&buffer, &int_val, TRUE);
	74	dem->utm_zone = int_val;
	75	/* TODO -- southern or northern hemisphere?! */
	76	dem->utm_letter = 'N';
	77
	78	/* skip numbers 5-19 */
	79	for ( i = 0; i < 15; i++ ) {
	80	if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
	81	g_warning ("Invalid DEM header");
	82	return FALSE;
	83	}
	84	}
	85
	86	/* number 20 -- horizontal unit code (utm/ll) */
	87	get_double_and_continue(&buffer, &val, TRUE);
	88	dem->horiz_units = val;
	89	get_double_and_continue(&buffer, &val, TRUE);
	90	/* dem->orig_vert_units = val; now done below */
	91
	92	/* TODO: do this for real. these are only for 1:24k and 1:250k USGS */
	93	if ( dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS ) {
	94	dem->east_scale = 10.0; /* meters */
	95	dem->north_scale = 10.0;
	96	dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
	97	} else {
	98	dem->east_scale = 3.0; /* arcseconds */
	99	dem->north_scale = 3.0;
	100	dem->orig_vert_units = VIK_DEM_VERT_METERS;
	101	}
	102
	103	/* skip next */
	104	get_double_and_continue(&buffer, &val, TRUE);
	105
	106	/* now we get the four corner points. record the min and max. */
	107	get_double_and_continue(&buffer, &val, TRUE);
	108	dem->min_east = dem->max_east = val;
	109	get_double_and_continue(&buffer, &val, TRUE);
	110	dem->min_north = dem->max_north = val;
	111
	112	for ( i = 0; i < 3; i++ ) {
	113	get_double_and_continue(&buffer, &val, TRUE);
	114	if ( val < dem->min_east ) dem->min_east = val;
	115	if ( val > dem->max_east ) dem->max_east = val;
	116	get_double_and_continue(&buffer, &val, TRUE);
	117	if ( val < dem->min_north ) dem->min_north = val;
	118	if ( val > dem->max_north ) dem->max_north = val;
	119	}
	120
	121	return TRUE;
	122	}
	123
	124	static void dem_parse_block_as_cont ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	125	{
	126	gint tmp;
	127	while ( cur_row < GET_COLUMN(dem, cur_column)->n_points ) {
	128	if ( get_int_and_continue(&buffer, &tmp,FALSE) ) {
	129	if ( dem->orig_vert_units == VIK_DEM_VERT_DECIMETERS )
	130	GET_COLUMN(dem, cur_column)->points[cur_row] = (gint16) (tmp / 10);
	131	else
	132	GET_COLUMN(dem, cur_column)->points[cur_row] = (gint16) tmp;
	133	} else
	134	return;
	135	(*cur_row)++;
	136	}
	137	cur_row = -1; / expecting new column */
	138	}
	139
	140	static void dem_parse_block_as_header ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	141	{
	142	guint n_rows;
	143	gint i;
	144	gdouble east_west, south;
	145	gdouble tmp;
	146
	147	/* 1 x n_rows 1 east_west south x x x DATA */
	148
	149	if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) \|\| tmp != 1 ) {
	150	g_warning("Incorrect DEM Class B record: expected 1");
	151	return;
	152	}
	153
	154	/* don't need this */
	155	if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;
	156
	157	/* n_rows */
	158	if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
	159	return;
	160	n_rows = (guint) tmp;
	161
	162	if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) \|\| tmp != 1 ) {
	163	g_warning("Incorrect DEM Class B record: expected 1");
	164	return;
	165	}
	166
	167	if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
	168	return;
	169	if ( !get_double_and_continue(&buffer, &south, TRUE) )
	170	return;
	171
	172	/* next three things we don't need */
	173	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	174	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	175	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	176
	177
	178	dem->n_columns ++;
	179	(*cur_column) ++;
	180
	181	/* empty spaces for things before that were skipped */
	182	(*cur_row) = (south - dem->min_north) / dem->north_scale;
	183	if ( south > dem->max_north \|\| (*cur_row) < 0 )
	184	(*cur_row) = 0;
	185
	186	n_rows += *cur_row;
	187
	188	g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
	189	GET_COLUMN(dem,*cur_column)->east_west = east_west;
	190	GET_COLUMN(dem,*cur_column)->south = south;
	191	GET_COLUMN(dem,*cur_column)->n_points = n_rows;
	192	GET_COLUMN(dem,cur_column)->points = g_malloc(sizeof(gint16)n_rows);
	193
	194	/* no information for things before that */
	195	for ( i = 0; i < (*cur_row); i++ )
	196	GET_COLUMN(dem,*cur_column)->points[i] = VIK_DEM_INVALID_ELEVATION;
	197
	198	/* now just continue */
	199	dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );
	200
	201
	202	}
	203
	204	static void dem_parse_block ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	205	{
	206	/* if haven't read anything or have read all items in a columns and are expecting a new column */
	207	if ( cur_column == -1 \|\| cur_row == -1 ) {
	208	dem_parse_block_as_header(buffer, dem, cur_column, cur_row);
	209	} else {
	210	dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
	211	}
	212	}
	213
	214	/* return size of unzip data or 0 if failed */
	215	/* can be made generic to uncompress zip, gzip, bzip2 data */
	216	static guint uncompress_data(void uncompressed_buffer, guint uncompressed_size, void compressed_data, guint compressed_size)
	217	{
	218	z_stream stream;
	219	int err;
	220
	221	stream.next_in = compressed_data;
	222	stream.avail_in = compressed_size;
	223	stream.next_out = uncompressed_buffer;
	224	stream.avail_out = uncompressed_size;
	225	stream.zalloc = (alloc_func)0;
	226	stream.zfree = (free_func)0;
	227	stream.opaque = (voidpf)0;
	228
	229	/* negative windowBits to inflateInit2 means "no header" */
	230	if ((err = inflateInit2(&stream, -MAX_WBITS)) != Z_OK) {
	231	g_warning("%s(): inflateInit2 failed", __PRETTY_FUNCTION__);
	232	return 0;
	233	}
	234
	235	err = inflate(&stream, Z_FINISH);
	236	if ((err != Z_OK) && (err != Z_STREAM_END) && stream.msg) {
	237	g_warning("%s() inflate failed err=%d \"%s\"", __PRETTY_FUNCTION__, err, stream.msg == NULL ? "unknown" : stream.msg);
	238	inflateEnd(&stream);
	239	return 0;
	240	}
	241
	242	inflateEnd(&stream);
	243	return(stream.total_out);
	244	}
	245
	246	static void unzip_hgt_file(gchar zip_file, gulong *unzip_size)
	247	{
	248	void *unzip_data = NULL;
	249	gchar *zip_data;
	250	struct _lfh {
	251	guint32 sig;
	252	guint16 extract_version;
	253	guint16 flags;
	254	guint16 comp_method;
	255	guint16 time;
	256	guint16 date;
	257	guint32 crc_32;
	258	guint32 compressed_size;
	259	guint32 uncompressed_size;
	260	guint16 filename_len;
	261	guint16 extra_field_len;
	262	} __attribute__ ((__packed__)) *local_file_header = NULL;
	263
	264
	265	local_file_header = (struct _lfh *) zip_file;
	266	if (local_file_header->sig != 0x04034b50) {
	267	g_warning("%s(): wrong format\n", __PRETTY_FUNCTION__);
	268	g_free(unzip_data);
	269	goto end;
	270	}
	271
	272	zip_data = zip_file + sizeof(struct _lfh) + local_file_header->filename_len + local_file_header->extra_field_len;
	273	unzip_data = g_malloc(local_file_header->uncompressed_size);
	274	gulong uncompressed_size = local_file_header->uncompressed_size;
	275
	276	if (!(*unzip_size = uncompress_data(unzip_data, uncompressed_size, zip_data, local_file_header->compressed_size))) {
	277	g_free(unzip_data);
	278	unzip_data = NULL;
	279	goto end;
	280	}
	281
	282	end:
	283	return(unzip_data);
	284	}
	285
	286	static VikDEM vik_dem_read_srtm_hgt(FILE f, const gchar *basename, gboolean zip)
	287	{
	288	gint i, j;
	289	VikDEM *dem;
	290	struct stat stat;
	291	off_t file_size;
	292	gint16 *dem_mem = NULL;
	293	gint16 *dem_file = NULL;
	294	const gint num_rows_3sec = 1201;
	295	const gint num_rows_1sec = 3601;
	296	gint num_rows;
	297	int fd = fileno(f);
	298	gint arcsec;
	299
	300	dem = g_malloc(sizeof(VikDEM));
	301
	302	dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
	303	dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
	304
	305	/* TODO */
	306	dem->min_north = atoi(basename+1) * 3600;
	307	dem->min_east = atoi(basename+4) * 3600;
	308	if ( basename[0] == 'S' )
	309	dem->min_north = - dem->min_north;
	310	if ( basename[3] == 'W' )
	311	dem->min_east = - dem->min_east;
	312
	313	dem->max_north = 3600 + dem->min_north;
	314	dem->max_east = 3600 + dem->min_east;
	315
	316	dem->columns = g_ptr_array_new();
	317	dem->n_columns = 0;
	318
	319	if (fstat(fd, &stat) == -1)
	320	g_error("%s(): fstat failed on %s\n", __PRETTY_FUNCTION__, basename);
	321	if ((dem_file = mmap(NULL, stat.st_size, PROT_READ, MAP_PRIVATE, fd, 0)) == (void *) -1)
	322	g_error("%s(): mmap failed on %s\n", __PRETTY_FUNCTION__, basename);
	323
	324	file_size = stat.st_size;
	325	dem_mem = dem_file;
	326	if (zip) {
	327	void *unzip_mem = NULL;
	328	gulong ucsize;
	329
	330	if ((unzip_mem = unzip_hgt_file((gchar *)dem_file, &ucsize)) == NULL) {
	331	munmap(dem_file, file_size);
	332	g_ptr_array_free(dem->columns, TRUE);
	333	g_free(dem);
	334	return NULL;
	335	}
	336
	337	dem_mem = unzip_mem;
	338	file_size = ucsize;
	339	}
	340
	341	if (file_size == (num_rows_3sec * num_rows_3sec * sizeof(gint16)))
	342	arcsec = 3;
	343	else if (file_size == (num_rows_1sec * num_rows_1sec * sizeof(gint16)))
	344	arcsec = 1;
	345	else {
	346	g_warning("%s(): file %s does not have right size\n", __PRETTY_FUNCTION__, basename);
	347	munmap(dem_file, file_size);
	348	g_ptr_array_free(dem->columns, TRUE);
	349	g_free(dem);
	350	return NULL;
	351	}
	352
	353	num_rows = (arcsec == 3) ? num_rows_3sec : num_rows_1sec;
	354	dem->east_scale = dem->north_scale = arcsec;
	355
	356	for ( i = 0; i < num_rows; i++ ) {
	357	dem->n_columns++;
	358	g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
	359	GET_COLUMN(dem,i)->east_west = dem->min_east + arcsec*i;
	360	GET_COLUMN(dem,i)->south = dem->min_north;
	361	GET_COLUMN(dem,i)->n_points = num_rows;
	362	GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*num_rows);
	363	}
	364
	365	int ent = 0;
	366	for ( i = (num_rows - 1); i >= 0; i-- ) {
	367	for ( j = 0; j < num_rows; j++ ) {
	368	GET_COLUMN(dem,j)->points[i] = GINT16_FROM_BE(dem_mem[ent]);
	369	ent++;
	370	}
	371
	372	}
	373
	374	if (zip)
	375	g_free(dem_mem);
	376	munmap(dem_file, stat.st_size);
	377	return dem;
	378	}
	379
	380	#define IS_SRTM_HGT(fn) (strlen((fn))==11 && (fn)[7]=='.' && (fn)[8]=='h' && (fn)[9]=='g' && (fn)[10]=='t' && ((fn)[0]=='N' \|\| (fn)[0]=='S') && ((fn)[3]=='E' \|\| (fn)[3]=='W'))
	381
	382	VikDEM vik_dem_new_from_file(const gchar file)
	383	{
	384	FILE *f;
	385	VikDEM *rv;
	386	gchar buffer[DEM_BLOCK_SIZE+1];
	387
	388	/* use to record state for dem_parse_block */
	389	gint cur_column = -1;
	390	gint cur_row = -1;
	391	const gchar *basename = a_file_basename(file);
	392
	393	/* FILE IO */
	394	f = fopen(file, "r");
	395	if ( !f )
	396	return NULL;
	397
	398	if ( (strlen(basename)==11 \|\| ((strlen(basename) == 15) && (basename[11] == '.' && basename[12] == 'z' && basename[13] == 'i' && basename[14] == 'p'))) &&
	399	basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
	400	(basename[0] == 'N' \|\| basename[0] == 'S') && (basename[3] == 'E' \|\| basename[3] =='W')) {
	401	gboolean is_zip_file = (strlen(basename) == 15);
	402	rv = vik_dem_read_srtm_hgt(f, basename, is_zip_file);
	403	fclose(f);
	404	return(rv);
	405	}
	406
	407	/* Create Structure */
	408	rv = g_malloc(sizeof(VikDEM));
	409
	410	/* Header */
	411	buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
	412	if ( ! dem_parse_header ( buffer, rv ) ) {
	413	g_free ( rv );
	414	return NULL;
	415	}
	416	/* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */
	417
	418	rv->columns = g_ptr_array_new();
	419	rv->n_columns = 0;
	420
	421	/* Column -- Data */
	422	while (! feof(f) ) {
	423	gchar *tmp;
	424
	425	/* read block */
	426	buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
	427
	428	/* Fix Fortran-style exponentiation */
	429	tmp = buffer;
	430	while (*tmp) {
	431	if ( *tmp=='D')
	432	*tmp='E';
	433	tmp++;
	434	}
	435
	436	dem_parse_block(buffer, rv, &cur_column, &cur_row);
	437	}
	438
	439	/* TODO - class C records (right now says 'Invalid' and dies) */
	440
	441	fclose(f);
	442
	443	/* 24k scale */
	444	if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
	445	rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;
	446
	447	/* FIXME bug in 10m DEM's */
	448	if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
	449	rv->min_east -= 100;
	450	rv->min_north += 200;
	451	}
	452
	453
	454	return rv;
	455	}
	456
	457	void vik_dem_free ( VikDEM *dem )
	458	{
	459	guint i;
	460	for ( i = 0; i < dem->n_columns; i++)
	461	g_free ( GET_COLUMN(dem, i)->points );
	462	g_ptr_array_free ( dem->columns, TRUE );
	463	g_free ( dem );
	464	}
	465
	466	gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
	467	{
	468	if ( col < dem->n_columns )
	469	if ( row < GET_COLUMN(dem, col)->n_points )
	470	return GET_COLUMN(dem, col)->points[row];
	471	return VIK_DEM_INVALID_ELEVATION;
	472	}
	473
	474	gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
	475	{
	476	gint col, row;
	477
	478	if ( east > dem->max_east \|\| east < dem->min_east \|\|
	479	north > dem->max_north \|\| north < dem->min_north )
	480	return VIK_DEM_INVALID_ELEVATION;
	481
	482	col = (gint) floor((east - dem->min_east) / dem->east_scale);
	483	row = (gint) floor((north - dem->min_north) / dem->north_scale);
	484
	485	return vik_dem_get_xy ( dem, col, row );
	486	}
	487
	488	static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
	489	gdouble east, gdouble north, /* in seconds */
	490	gint16 elevs, gint16 dists)
	491	{
	492	int i;
	493	int cols[4], rows[4];
	494	struct LatLon ll[4];
	495	struct LatLon pos;
	496
	497	if ( east > dem->max_east \|\| east < dem->min_east \|\|
	498	north > dem->max_north \|\| north < dem->min_north )
	499	return FALSE; /* got nothing */
	500
	501	pos.lon = east/3600;
	502	pos.lat = north/3600;
	503
	504	/* order of the data: sw, nw, ne, se */
	505	cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale); /* sw */
	506	rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
	507	ll[0].lon = (east - fabs(fmod(east, dem->east_scale)))/3600;
	508	ll[0].lat = (north - fabs(fmod(north, dem->north_scale)))/3600;
	509
	510	cols[1] = cols[0]; /nw/
	511	rows[1] = rows[0] + 1;
	512	ll[1].lon = ll[0].lon;
	513	ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
	514
	515	cols[2] = cols[0] + 1; /ne/
	516	rows[2] = rows[0] + 1;
	517	ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
	518	ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
	519
	520	cols[3] = cols[0] + 1; /se/
	521	rows[3] = rows[0];
	522	ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
	523	ll[3].lat = ll[0].lat;
	524
	525	for (i = 0; i < 4; i++) {
	526	if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION)
	527	return FALSE;
	528	dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
	529	}
	530
	531	return TRUE; /* all OK */
	532	}
	533
	534	gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
	535	{
	536	int i;
	537	gint16 elevs[4], dists[4];
	538
	539	if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
	540	return VIK_DEM_INVALID_ELEVATION;
	541
	542	for (i = 0; i < 4; i++) {
	543	if (dists[i] < 1) {
	544	return(elevs[i]);
	545	}
	546	}
	547
	548	gdouble t = (gdouble)elevs[0]/dists[0] + (gdouble)elevs[1]/dists[1] + (gdouble)elevs[2]/dists[2] + (gdouble)elevs[3]/dists[3];
	549	gdouble b = 1.0/dists[0] + 1.0/dists[1] + 1.0/dists[2] + 1.0/dists[3];
	550
	551	return(t/b);
	552	}
	553
	554	gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
	555	{
	556	return vik_dem_get_simple_interpol(dem, east, north);
	557	}
	558
	559	void vik_dem_east_north_to_xy ( VikDEM dem, gdouble east, gdouble north, guint col, guint *row )
	560	{
	561	*col = (gint) floor((east - dem->min_east) / dem->east_scale);
	562	*row = (gint) floor((north - dem->min_north) / dem->north_scale);
	563	if ( col < 0 ) col = 0;
	564	if ( row < 0 ) row = 0;
	565	}
	566