[andy/viking.git] / src / dem.c

/*
 * viking -- GPS Data and Topo Analyzer, Explorer, and Manager
 *
 * Copyright (C) 2003-2008, Evan Battaglia <gtoevan@gmx.net>
 * Copyright (C) 2007, Quy Tonthat <qtonthat@gmail.com>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#ifdef HAVE_STRING_H
#include <string.h>
#endif
#include <glib.h>
#ifdef HAVE_MATH_H
#include <math.h>
#endif
#ifdef HAVE_STDLIB_H
#include <stdlib.h>
#endif
#include <zlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif

#include <glib/gstdio.h>
#include <glib/gi18n.h>

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

/* Compatibility */
#if ! GLIB_CHECK_VERSION(2,22,0)
#define g_mapped_file_unref g_mapped_file_free
#endif

#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 (GUINT32_FROM_LE(local_file_header->sig) != 0x04034b50) {
		g_warning("%s(): wrong format", __PRETTY_FUNCTION__);
		g_free(unzip_data);
		goto end;
	}

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

	if (!(*unzip_size = uncompress_data(unzip_data, uncompressed_size, zip_data, GUINT32_FROM_LE(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(const gchar *file_name, const gchar *basename, gboolean zip)
{
  gint i, j;
  VikDEM *dem;
  off_t file_size;
  gint16 *dem_mem = NULL;
  gchar *dem_file = NULL;
  const gint num_rows_3sec = 1201;
  const gint num_rows_1sec = 3601;
  gint num_rows;
  GMappedFile *mf;
  gint arcsec;
  GError *error = NULL;

  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 ((mf = g_mapped_file_new(file_name, FALSE, &error)) == NULL) {
    g_critical(_("Couldn't map file %s: %s"), file_name, error->message);
    g_error_free(error);
    g_free(dem);
    return NULL;
  }
  file_size = g_mapped_file_get_length(mf);
  dem_file = g_mapped_file_get_contents(mf);
  
  if (zip) {
    void *unzip_mem = NULL;
    gulong ucsize;

    if ((unzip_mem = unzip_hgt_file(dem_file, &ucsize)) == NULL) {
      g_mapped_file_unref(mf);
      g_ptr_array_free(dem->columns, TRUE);
      g_free(dem);
      return NULL;
    }

    dem_mem = unzip_mem;
    file_size = ucsize;
  }
  else
    dem_mem = (gint16 *)dem_file;

  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", __PRETTY_FUNCTION__, basename);
    g_mapped_file_unref(mf);
    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);
  g_mapped_file_unref(mf);
  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=NULL;
  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);

  if ( g_access ( file, R_OK ) != 0 )
    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(file, basename, is_zip_file);
    return(rv);
  }

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

      /* Header */
  f = g_fopen(file, "r");
  if ( !f ) {
    g_free ( rv );
    return NULL;
  }
  buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
  if ( ! dem_parse_header ( buffer, rv ) ) {
    g_free ( rv );
    fclose(f);
    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);
  f = NULL;

  /* 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 */
  /* sw */
  cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);
  rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
  ll[0].lon = (dem->min_east + dem->east_scale*cols[0])/3600;
  ll[0].lat = (dem->min_north + dem->north_scale*rows[0])/3600;
  /* nw */
  cols[1] = cols[0];
  rows[1] = rows[0] + 1;
  ll[1].lon = ll[0].lon;
  ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
  /* ne */
  cols[2] = cols[0] + 1;
  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;
  /* se */
  cols[3] = cols[0] + 1;
  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]);
  }

#if 0  /* debug */
  for (i = 0; i < 4; i++)
    fprintf(stderr, "%f:%f:%d:%d  ", ll[i].lat, ll[i].lon, dists[i], elevs[i]);
  fprintf(stderr, "   north_scale=%f\n", dem->north_scale);
#endif

  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 )
{
  int i;
  gint16 elevs[4], dists[4];
  gint16 max_dist;
  gdouble t = 0.0;
  gdouble b = 0.0;

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

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

  gdouble tmp;
#if 0 /* derived method by Franke & Nielson. Does not seem to work too well here */
  for (i = 0; i < 4; i++) {
    tmp = pow((1.0*(max_dist - dists[i])/max_dist*dists[i]), 2);
    t += tmp*elevs[i];
    b += tmp;
  }
#endif

  for (i = 0; i < 4; i++) {
    tmp = pow((1.0/dists[i]), 2);
    t += tmp*elevs[i];
    b += tmp;
  }

  // fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);

  return(t/b);

}

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
f6171fe3 EB	1	/*
	2	* viking -- GPS Data and Topo Analyzer, Explorer, and Manager
	3	*
	4	* Copyright (C) 2003-2008, Evan Battaglia <gtoevan@gmx.net>
a482007a	5	* Copyright (C) 2007, Quy Tonthat <qtonthat@gmail.com>
f6171fe3 EB	6	*
	7	* This program is free software; you can redistribute it and/or modify
	8	* it under the terms of the GNU General Public License as published by
	9	* the Free Software Foundation; either version 2 of the License, or
	10	* (at your option) any later version.
	11	*
	12	* This program is distributed in the hope that it will be useful,
	13	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	14	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	15	* GNU General Public License for more details.
	16	*
	17	* You should have received a copy of the GNU General Public License
	18	* along with this program; if not, write to the Free Software
	19	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
	20	*
	21	*/
	22
72b6a3ed JB	23	#ifdef HAVE_CONFIG_H
	24	#include "config.h"
	25	#endif
	26
ad0a8c2d	27	#include <stdio.h>
8c00358d	28	#ifdef HAVE_STRING_H
ad0a8c2d	29	#include <string.h>
8c00358d	30	#endif
ad0a8c2d	31	#include <glib.h>
8c00358d	32	#ifdef HAVE_MATH_H
ad0a8c2d	33	#include <math.h>
8c00358d GB	34	#endif
8c00358d GB	35	#ifdef HAVE_STDLIB_H
ad0a8c2d	36	#include <stdlib.h>
8c00358d	37	#endif
46e4fa54 QT	38	#include <zlib.h>
	39	#include <sys/types.h>
	40	#include <sys/stat.h>
31766324 MA	41	#ifdef HAVE_UNISTD_H
	42	#include <unistd.h>
	43	#endif
46e4fa54	44
6af463da	45	#include <glib/gstdio.h>
4c77d5e0 GB	46	#include <glib/gi18n.h>
4c77d5e0 GB	47
ad0a8c2d EB	48	#include "dem.h"
	49	#include "file.h"
	50
316945d8 GB	51	/* Compatibility */
	52	#if ! GLIB_CHECK_VERSION(2,22,0)
	53	#define g_mapped_file_unref g_mapped_file_free
	54	#endif
	55
ad0a8c2d EB	56	#define DEM_BLOCK_SIZE 1024
	57	#define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))
	58
	59	static gboolean get_double_and_continue ( gchar *buffer, gdouble tmp, gboolean warn )
	60	{
	61	gchar *endptr;
	62	tmp = g_strtod(buffer, &endptr);
	63	if ( endptr == NULL\|\| endptr == *buffer ) {
	64	if ( warn )
4c77d5e0	65	g_warning(_("Invalid DEM"));
ad0a8c2d EB	66	return FALSE;
	67	}
	68	*buffer=endptr;
	69	return TRUE;
	70	}
	71
	72
	73	static gboolean get_int_and_continue ( gchar *buffer, gint tmp, gboolean warn )
	74	{
	75	gchar *endptr;
	76	tmp = strtol(buffer, &endptr, 10);
	77	if ( endptr == NULL\|\| endptr == *buffer ) {
	78	if ( warn )
4c77d5e0	79	g_warning(_("Invalid DEM"));
ad0a8c2d EB	80	return FALSE;
	81	}
	82	*buffer=endptr;
	83	return TRUE;
	84	}
	85
	86	static gboolean dem_parse_header ( gchar buffer, VikDEM dem )
	87	{
	88	gdouble val;
	89	gint int_val;
	90	guint i;
	91	gchar *tmp = buffer;
	92
	93	/* incomplete header */
	94	if ( strlen(buffer) != 1024 )
	95	return FALSE;
	96
	97	/* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
	98	while (*tmp) {
	99	if ( *tmp=='D')
	100	*tmp='E';
	101	tmp++;
	102	}
	103
	104	/* skip name */
	105	buffer += 149;
	106
	107	/* "DEM level code, pattern code, palaimetric reference system code" -- skip */
	108	get_int_and_continue(&buffer, &int_val, TRUE);
	109	get_int_and_continue(&buffer, &int_val, TRUE);
	110	get_int_and_continue(&buffer, &int_val, TRUE);
	111
	112	/* zone */
	113	get_int_and_continue(&buffer, &int_val, TRUE);
	114	dem->utm_zone = int_val;
	115	/* TODO -- southern or northern hemisphere?! */
	116	dem->utm_letter = 'N';
	117
	118	/* skip numbers 5-19 */
	119	for ( i = 0; i < 15; i++ ) {
	120	if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
4c77d5e0	121	g_warning (_("Invalid DEM header"));
ad0a8c2d EB	122	return FALSE;
	123	}
	124	}
	125
	126	/* number 20 -- horizontal unit code (utm/ll) */
	127	get_double_and_continue(&buffer, &val, TRUE);
	128	dem->horiz_units = val;
	129	get_double_and_continue(&buffer, &val, TRUE);
70147e60	130	/* dem->orig_vert_units = val; now done below */
ad0a8c2d EB	131
	132	/* TODO: do this for real. these are only for 1:24k and 1:250k USGS */
	133	if ( dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS ) {
	134	dem->east_scale = 10.0; /* meters */
	135	dem->north_scale = 10.0;
70147e60	136	dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
ad0a8c2d EB	137	} else {
	138	dem->east_scale = 3.0; /* arcseconds */
	139	dem->north_scale = 3.0;
70147e60	140	dem->orig_vert_units = VIK_DEM_VERT_METERS;
ad0a8c2d EB	141	}
	142
	143	/* skip next */
	144	get_double_and_continue(&buffer, &val, TRUE);
	145
	146	/* now we get the four corner points. record the min and max. */
	147	get_double_and_continue(&buffer, &val, TRUE);
	148	dem->min_east = dem->max_east = val;
	149	get_double_and_continue(&buffer, &val, TRUE);
	150	dem->min_north = dem->max_north = val;
	151
	152	for ( i = 0; i < 3; i++ ) {
	153	get_double_and_continue(&buffer, &val, TRUE);
	154	if ( val < dem->min_east ) dem->min_east = val;
	155	if ( val > dem->max_east ) dem->max_east = val;
	156	get_double_and_continue(&buffer, &val, TRUE);
	157	if ( val < dem->min_north ) dem->min_north = val;
	158	if ( val > dem->max_north ) dem->max_north = val;
	159	}
	160
	161	return TRUE;
	162	}
	163
	164	static void dem_parse_block_as_cont ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	165	{
	166	gint tmp;
	167	while ( cur_row < GET_COLUMN(dem, cur_column)->n_points ) {
70147e60 EB	168	if ( get_int_and_continue(&buffer, &tmp,FALSE) ) {
	169	if ( dem->orig_vert_units == VIK_DEM_VERT_DECIMETERS )
	170	GET_COLUMN(dem, cur_column)->points[cur_row] = (gint16) (tmp / 10);
	171	else
	172	GET_COLUMN(dem, cur_column)->points[cur_row] = (gint16) tmp;
	173	} else
ad0a8c2d EB	174	return;
	175	(*cur_row)++;
	176	}
	177	cur_row = -1; / expecting new column */
	178	}
	179
	180	static void dem_parse_block_as_header ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	181	{
	182	guint n_rows;
	183	gint i;
	184	gdouble east_west, south;
	185	gdouble tmp;
	186
	187	/* 1 x n_rows 1 east_west south x x x DATA */
	188
	189	if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) \|\| tmp != 1 ) {
4c77d5e0	190	g_warning(_("Incorrect DEM Class B record: expected 1"));
ad0a8c2d EB	191	return;
	192	}
	193
	194	/* don't need this */
	195	if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;
	196
	197	/* n_rows */
	198	if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
	199	return;
	200	n_rows = (guint) tmp;
	201
	202	if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) \|\| tmp != 1 ) {
4c77d5e0	203	g_warning(_("Incorrect DEM Class B record: expected 1"));
ad0a8c2d EB	204	return;
	205	}
	206
	207	if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
	208	return;
	209	if ( !get_double_and_continue(&buffer, &south, TRUE) )
	210	return;
	211
	212	/* next three things we don't need */
	213	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	214	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	215	if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
	216
	217
	218	dem->n_columns ++;
	219	(*cur_column) ++;
	220
	221	/* empty spaces for things before that were skipped */
	222	(*cur_row) = (south - dem->min_north) / dem->north_scale;
	223	if ( south > dem->max_north \|\| (*cur_row) < 0 )
	224	(*cur_row) = 0;
	225
	226	n_rows += *cur_row;
	227
	228	g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
	229	GET_COLUMN(dem,*cur_column)->east_west = east_west;
	230	GET_COLUMN(dem,*cur_column)->south = south;
	231	GET_COLUMN(dem,*cur_column)->n_points = n_rows;
	232	GET_COLUMN(dem,cur_column)->points = g_malloc(sizeof(gint16)n_rows);
	233
	234	/* no information for things before that */
	235	for ( i = 0; i < (*cur_row); i++ )
	236	GET_COLUMN(dem,*cur_column)->points[i] = VIK_DEM_INVALID_ELEVATION;
	237
	238	/* now just continue */
	239	dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );
	240
	241
	242	}
	243
	244	static void dem_parse_block ( gchar buffer, VikDEM dem, gint cur_column, gint cur_row )
	245	{
	246	/* if haven't read anything or have read all items in a columns and are expecting a new column */
	247	if ( cur_column == -1 \|\| cur_row == -1 ) {
	248	dem_parse_block_as_header(buffer, dem, cur_column, cur_row);
	249	} else {
	250	dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
	251	}
	252	}
	253
46e4fa54 QT	254	/* return size of unzip data or 0 if failed */
	255	/* can be made generic to uncompress zip, gzip, bzip2 data */
	256	static guint uncompress_data(void uncompressed_buffer, guint uncompressed_size, void compressed_data, guint compressed_size)
	257	{
	258	z_stream stream;
	259	int err;
	260
	261	stream.next_in = compressed_data;
	262	stream.avail_in = compressed_size;
	263	stream.next_out = uncompressed_buffer;
	264	stream.avail_out = uncompressed_size;
	265	stream.zalloc = (alloc_func)0;
	266	stream.zfree = (free_func)0;
	267	stream.opaque = (voidpf)0;
	268
	269	/* negative windowBits to inflateInit2 means "no header" */
	270	if ((err = inflateInit2(&stream, -MAX_WBITS)) != Z_OK) {
	271	g_warning("%s(): inflateInit2 failed", __PRETTY_FUNCTION__);
	272	return 0;
	273	}
	274
	275	err = inflate(&stream, Z_FINISH);
	276	if ((err != Z_OK) && (err != Z_STREAM_END) && stream.msg) {
	277	g_warning("%s() inflate failed err=%d \"%s\"", __PRETTY_FUNCTION__, err, stream.msg == NULL ? "unknown" : stream.msg);
	278	inflateEnd(&stream);
	279	return 0;
	280	}
	281
	282	inflateEnd(&stream);
	283	return(stream.total_out);
	284	}
	285
	286	static void unzip_hgt_file(gchar zip_file, gulong *unzip_size)
	287	{
	288	void *unzip_data = NULL;
	289	gchar *zip_data;
	290	struct _lfh {
	291	guint32 sig;
	292	guint16 extract_version;
	293	guint16 flags;
	294	guint16 comp_method;
	295	guint16 time;
	296	guint16 date;
	297	guint32 crc_32;
	298	guint32 compressed_size;
	299	guint32 uncompressed_size;
	300	guint16 filename_len;
	301	guint16 extra_field_len;
	302	} __attribute__ ((__packed__)) *local_file_header = NULL;
	303
	304
	305	local_file_header = (struct _lfh *) zip_file;
0a8cbd7b	306	if (GUINT32_FROM_LE(local_file_header->sig) != 0x04034b50) {
4258f4e2	307	g_warning("%s(): wrong format", __PRETTY_FUNCTION__);
46e4fa54 QT	308	g_free(unzip_data);
	309	goto end;
	310	}
	311
0a8cbd7b MK	312	zip_data = zip_file + sizeof(struct _lfh)
	313	+ GUINT16_FROM_LE(local_file_header->filename_len)
	314	+ GUINT16_FROM_LE(local_file_header->extra_field_len);
	315	gulong uncompressed_size = GUINT32_FROM_LE(local_file_header->uncompressed_size);
	316	unzip_data = g_malloc(uncompressed_size);
46e4fa54	317
0a8cbd7b	318	if (!(*unzip_size = uncompress_data(unzip_data, uncompressed_size, zip_data, GUINT32_FROM_LE(local_file_header->compressed_size)))) {
46e4fa54 QT	319	g_free(unzip_data);
	320	unzip_data = NULL;
	321	goto end;
	322	}
	323
	324	end:
	325	return(unzip_data);
	326	}
	327
31766324	328	static VikDEM vik_dem_read_srtm_hgt(const gchar file_name, const gchar *basename, gboolean zip)
ad0a8c2d	329	{
ad0a8c2d EB	330	gint i, j;
ad0a8c2d EB	331	VikDEM *dem;
46e4fa54 QT	332	off_t file_size;
46e4fa54 QT	333	gint16 *dem_mem = NULL;
31766324	334	gchar *dem_file = NULL;
46e4fa54 QT	335	const gint num_rows_3sec = 1201;
	336	const gint num_rows_1sec = 3601;
	337	gint num_rows;
31766324	338	GMappedFile *mf;
46e4fa54	339	gint arcsec;
a95be497	340	GError *error = NULL;
ad0a8c2d EB	341
	342	dem = g_malloc(sizeof(VikDEM));
	343
	344	dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
70147e60	345	dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
ad0a8c2d EB	346
	347	/* TODO */
	348	dem->min_north = atoi(basename+1) * 3600;
	349	dem->min_east = atoi(basename+4) * 3600;
	350	if ( basename[0] == 'S' )
	351	dem->min_north = - dem->min_north;
	352	if ( basename[3] == 'W' )
	353	dem->min_east = - dem->min_east;
	354
	355	dem->max_north = 3600 + dem->min_north;
	356	dem->max_east = 3600 + dem->min_east;
	357
	358	dem->columns = g_ptr_array_new();
	359	dem->n_columns = 0;
	360
a95be497	361	if ((mf = g_mapped_file_new(file_name, FALSE, &error)) == NULL) {
24ff9c7b	362	g_critical(_("Couldn't map file %s: %s"), file_name, error->message);
a95be497	363	g_error_free(error);
31766324 MA	364	g_free(dem);
	365	return NULL;
	366	}
	367	file_size = g_mapped_file_get_length(mf);
	368	dem_file = g_mapped_file_get_contents(mf);
	369
46e4fa54 QT	370	if (zip) {
	371	void *unzip_mem = NULL;
	372	gulong ucsize;
	373
31766324	374	if ((unzip_mem = unzip_hgt_file(dem_file, &ucsize)) == NULL) {
316945d8	375	g_mapped_file_unref(mf);
46e4fa54 QT	376	g_ptr_array_free(dem->columns, TRUE);
	377	g_free(dem);
	378	return NULL;
	379	}
	380
	381	dem_mem = unzip_mem;
	382	file_size = ucsize;
	383	}
92390c64 CK	384	else
92390c64 CK	385	dem_mem = (gint16 *)dem_file;
46e4fa54 QT	386
	387	if (file_size == (num_rows_3sec * num_rows_3sec * sizeof(gint16)))
	388	arcsec = 3;
	389	else if (file_size == (num_rows_1sec * num_rows_1sec * sizeof(gint16)))
	390	arcsec = 1;
	391	else {
4258f4e2	392	g_warning("%s(): file %s does not have right size", __PRETTY_FUNCTION__, basename);
316945d8	393	g_mapped_file_unref(mf);
46e4fa54 QT	394	g_free(dem);
	395	return NULL;
	396	}
	397
	398	num_rows = (arcsec == 3) ? num_rows_3sec : num_rows_1sec;
	399	dem->east_scale = dem->north_scale = arcsec;
	400
	401	for ( i = 0; i < num_rows; i++ ) {
ad0a8c2d EB	402	dem->n_columns++;
ad0a8c2d EB	403	g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
46e4fa54	404	GET_COLUMN(dem,i)->east_west = dem->min_east + arcsec*i;
ad0a8c2d	405	GET_COLUMN(dem,i)->south = dem->min_north;
46e4fa54 QT	406	GET_COLUMN(dem,i)->n_points = num_rows;
46e4fa54 QT	407	GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*num_rows);
ad0a8c2d EB	408	}
ad0a8c2d EB	409
46e4fa54 QT	410	int ent = 0;
	411	for ( i = (num_rows - 1); i >= 0; i-- ) {
	412	for ( j = 0; j < num_rows; j++ ) {
792a7488 QT	413	GET_COLUMN(dem,j)->points[i] = GINT16_FROM_BE(dem_mem[ent]);
792a7488 QT	414	ent++;
ad0a8c2d EB	415	}
	416
	417	}
	418
46e4fa54 QT	419	if (zip)
46e4fa54 QT	420	g_free(dem_mem);
316945d8	421	g_mapped_file_unref(mf);
ad0a8c2d EB	422	return dem;
	423	}
	424
	425	#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'))
	426
	427	VikDEM vik_dem_new_from_file(const gchar file)
	428	{
31766324	429	FILE *f=NULL;
ad0a8c2d EB	430	VikDEM *rv;
	431	gchar buffer[DEM_BLOCK_SIZE+1];
	432
	433	/* use to record state for dem_parse_block */
	434	gint cur_column = -1;
	435	gint cur_row = -1;
	436	const gchar *basename = a_file_basename(file);
	437
31766324	438	if ( g_access ( file, R_OK ) != 0 )
ad0a8c2d EB	439	return NULL;
ad0a8c2d EB	440
46e4fa54 QT	441	if ( (strlen(basename)==11 \|\| ((strlen(basename) == 15) && (basename[11] == '.' && basename[12] == 'z' && basename[13] == 'i' && basename[14] == 'p'))) &&
46e4fa54 QT	442	basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
2dc9ab5e	443	(basename[0] == 'N' \|\| basename[0] == 'S') && (basename[3] == 'E' \|\| basename[3] =='W')) {
46e4fa54	444	gboolean is_zip_file = (strlen(basename) == 15);
31766324	445	rv = vik_dem_read_srtm_hgt(file, basename, is_zip_file);
2dc9ab5e QT	446	return(rv);
2dc9ab5e QT	447	}
ad0a8c2d EB	448
	449	/* Create Structure */
	450	rv = g_malloc(sizeof(VikDEM));
	451
	452	/* Header */
31766324	453	f = g_fopen(file, "r");
3b36279c JJ	454	if ( !f ) {
3b36279c JJ	455	g_free ( rv );
31766324	456	return NULL;
3b36279c	457	}
ad0a8c2d EB	458	buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
	459	if ( ! dem_parse_header ( buffer, rv ) ) {
	460	g_free ( rv );
dad45413	461	fclose(f);
ad0a8c2d EB	462	return NULL;
	463	}
	464	/* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */
	465
	466	rv->columns = g_ptr_array_new();
	467	rv->n_columns = 0;
	468
	469	/* Column -- Data */
	470	while (! feof(f) ) {
	471	gchar *tmp;
	472
	473	/* read block */
	474	buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
	475
	476	/* Fix Fortran-style exponentiation */
	477	tmp = buffer;
	478	while (*tmp) {
	479	if ( *tmp=='D')
	480	*tmp='E';
	481	tmp++;
	482	}
	483
	484	dem_parse_block(buffer, rv, &cur_column, &cur_row);
	485	}
	486
	487	/* TODO - class C records (right now says 'Invalid' and dies) */
	488
	489	fclose(f);
8c060406	490	f = NULL;
ad0a8c2d EB	491
	492	/* 24k scale */
	493	if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
	494	rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;
	495
	496	/* FIXME bug in 10m DEM's */
	497	if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
	498	rv->min_east -= 100;
	499	rv->min_north += 200;
	500	}
	501
	502
	503	return rv;
	504	}
	505
	506	void vik_dem_free ( VikDEM *dem )
	507	{
	508	guint i;
	509	for ( i = 0; i < dem->n_columns; i++)
	510	g_free ( GET_COLUMN(dem, i)->points );
	511	g_ptr_array_free ( dem->columns, TRUE );
	512	g_free ( dem );
	513	}
	514
	515	gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
	516	{
	517	if ( col < dem->n_columns )
	518	if ( row < GET_COLUMN(dem, col)->n_points )
	519	return GET_COLUMN(dem, col)->points[row];
	520	return VIK_DEM_INVALID_ELEVATION;
	521	}
	522
	523	gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
	524	{
	525	gint col, row;
	526
	527	if ( east > dem->max_east \|\| east < dem->min_east \|\|
	528	north > dem->max_north \|\| north < dem->min_north )
	529	return VIK_DEM_INVALID_ELEVATION;
	530
	531	col = (gint) floor((east - dem->min_east) / dem->east_scale);
	532	row = (gint) floor((north - dem->min_north) / dem->north_scale);
	533
	534	return vik_dem_get_xy ( dem, col, row );
	535	}
	536
5ef1d57e QT	537	static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
	538	gdouble east, gdouble north, /* in seconds */
	539	gint16 elevs, gint16 dists)
	540	{
	541	int i;
	542	int cols[4], rows[4];
	543	struct LatLon ll[4];
	544	struct LatLon pos;
	545
	546	if ( east > dem->max_east \|\| east < dem->min_east \|\|
	547	north > dem->max_north \|\| north < dem->min_north )
	548	return FALSE; /* got nothing */
	549
	550	pos.lon = east/3600;
	551	pos.lat = north/3600;
	552
	553	/* order of the data: sw, nw, ne, se */
245575e1 QT	554	/* sw */
245575e1 QT	555	cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);
5ef1d57e	556	rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
228213c5 QT	557	ll[0].lon = (dem->min_east + dem->east_scale*cols[0])/3600;
228213c5 QT	558	ll[0].lat = (dem->min_north + dem->north_scale*rows[0])/3600;
245575e1 QT	559	/* nw */
245575e1 QT	560	cols[1] = cols[0];
5ef1d57e QT	561	rows[1] = rows[0] + 1;
	562	ll[1].lon = ll[0].lon;
	563	ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
245575e1 QT	564	/* ne */
245575e1 QT	565	cols[2] = cols[0] + 1;
5ef1d57e QT	566	rows[2] = rows[0] + 1;
	567	ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
	568	ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
245575e1 QT	569	/* se */
245575e1 QT	570	cols[3] = cols[0] + 1;
5ef1d57e QT	571	rows[3] = rows[0];
	572	ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
	573	ll[3].lat = ll[0].lat;
	574
	575	for (i = 0; i < 4; i++) {
	576	if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION)
	577	return FALSE;
	578	dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
	579	}
	580
228213c5 QT	581	#if 0 /* debug */
	582	for (i = 0; i < 4; i++)
	583	fprintf(stderr, "%f:%f:%d:%d ", ll[i].lat, ll[i].lon, dists[i], elevs[i]);
	584	fprintf(stderr, " north_scale=%f\n", dem->north_scale);
	585	#endif
	586
5ef1d57e QT	587	return TRUE; /* all OK */
	588	}
	589
	590	gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
	591	{
	592	int i;
	593	gint16 elevs[4], dists[4];
	594
	595	if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
	596	return VIK_DEM_INVALID_ELEVATION;
	597
	598	for (i = 0; i < 4; i++) {
	599	if (dists[i] < 1) {
	600	return(elevs[i]);
	601	}
	602	}
	603
	604	gdouble t = (gdouble)elevs[0]/dists[0] + (gdouble)elevs[1]/dists[1] + (gdouble)elevs[2]/dists[2] + (gdouble)elevs[3]/dists[3];
	605	gdouble b = 1.0/dists[0] + 1.0/dists[1] + 1.0/dists[2] + 1.0/dists[3];
	606
	607	return(t/b);
	608	}
	609
	610	gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
	611	{
228213c5 QT	612	int i;
	613	gint16 elevs[4], dists[4];
	614	gint16 max_dist;
	615	gdouble t = 0.0;
	616	gdouble b = 0.0;
	617
	618	if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
	619	return VIK_DEM_INVALID_ELEVATION;
	620
	621	max_dist = 0;
	622	for (i = 0; i < 4; i++) {
	623	if (dists[i] < 1) {
	624	return(elevs[i]);
	625	}
	626	if (dists[i] > max_dist)
	627	max_dist = dists[i];
	628	}
	629
	630	gdouble tmp;
	631	#if 0 /* derived method by Franke & Nielson. Does not seem to work too well here */
	632	for (i = 0; i < 4; i++) {
	633	tmp = pow((1.0(max_dist - dists[i])/max_distdists[i]), 2);
	634	t += tmp*elevs[i];
	635	b += tmp;
	636	}
	637	#endif
	638
	639	for (i = 0; i < 4; i++) {
	640	tmp = pow((1.0/dists[i]), 2);
	641	t += tmp*elevs[i];
	642	b += tmp;
	643	}
	644
	645	// fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);
	646
	647	return(t/b);
	648
5ef1d57e QT	649	}
5ef1d57e QT	650
ad0a8c2d EB	651	void vik_dem_east_north_to_xy ( VikDEM dem, gdouble east, gdouble north, guint col, guint *row )
	652	{
	653	*col = (gint) floor((east - dem->min_east) / dem->east_scale);
	654	*row = (gint) floor((north - dem->min_north) / dem->north_scale);
	655	if ( col < 0 ) col = 0;
	656	if ( row < 0 ) row = 0;
	657	}
	658