src/dem.c

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