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  * Copyright (C) 2007, Quy Tonthat <qtonthat@gmail.com>
   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
  23 #ifdef HAVE_CONFIG_H
  24 #include "config.h"
  25 #endif
  26
  27 #include <stdio.h>
  28 #ifdef HAVE_STRING_H
  29 #include <string.h>
  30 #endif
  31 #include <glib.h>
  32 #ifdef HAVE_MATH_H
  33 #include <math.h>
  34 #endif
  35 #ifdef HAVE_STDLIB_H
  36 #include <stdlib.h>
  37 #endif
  38 #include <zlib.h>
  39 #include <sys/types.h>
  40 #include <sys/stat.h>
  41 #ifdef HAVE_UNISTD_H
  42 #include <unistd.h>
  43 #endif
  44
  45 #include <glib/gstdio.h>
  46 #include <glib/gi18n.h>
  47
  48 #include "dem.h"
  49 #include "file.h"
  50
  51 /* Compatibility */
  52 #if ! GLIB_CHECK_VERSION(2,22,0)
  53 #define g_mapped_file_unref g_mapped_file_free
  54 #endif
  55
  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 )
  65       g_warning(_("Invalid DEM"));
  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 )
  79       g_warning(_("Invalid DEM"));
  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) ) {
 121       g_warning (_("Invalid DEM header"));
 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);
 130   /* dem->orig_vert_units = val; now done below */
 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;
 136     dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
 137   } else {
 138     dem->east_scale = 3.0; /* arcseconds */
 139     dem->north_scale = 3.0;
 140     dem->orig_vert_units = VIK_DEM_VERT_METERS;
 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 ) {
 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
 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 ) {
 190     g_warning(_("Incorrect DEM Class B record: expected 1"));
 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 ) {
 203     g_warning(_("Incorrect DEM Class B record: expected 1"));
 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
 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;
 306         if (GUINT32_FROM_LE(local_file_header->sig) != 0x04034b50) {
 307                 g_warning("%s(): wrong format", __PRETTY_FUNCTION__);
 308                 g_free(unzip_data);
 309                 goto end;
 310         }
 311
 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);
 317
 318         if (!(*unzip_size = uncompress_data(unzip_data, uncompressed_size, zip_data, GUINT32_FROM_LE(local_file_header->compressed_size)))) {
 319                 g_free(unzip_data);
 320                 unzip_data = NULL;
 321                 goto end;
 322         }
 323
 324 end:
 325         return(unzip_data);
 326 }
 327
 328 static VikDEM *vik_dem_read_srtm_hgt(const gchar *file_name, const gchar *basename, gboolean zip)
 329 {
 330   gint i, j;
 331   VikDEM *dem;
 332   off_t file_size;
 333   gint16 *dem_mem = NULL;
 334   gchar *dem_file = NULL;
 335   const gint num_rows_3sec = 1201;
 336   const gint num_rows_1sec = 3601;
 337   gint num_rows;
 338   GMappedFile *mf;
 339   gint arcsec;
 340   GError *error = NULL;
 341
 342   dem = g_malloc(sizeof(VikDEM));
 343
 344   dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
 345   dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
 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
 361   if ((mf = g_mapped_file_new(file_name, FALSE, &error)) == NULL) {
 362     g_critical(_("Couldn't map file %s: %s"), file_name, error->message);
 363     g_error_free(error);
 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
 370   if (zip) {
 371     void *unzip_mem = NULL;
 372     gulong ucsize;
 373
 374     if ((unzip_mem = unzip_hgt_file(dem_file, &ucsize)) == NULL) {
 375       g_mapped_file_unref(mf);
 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   }
 384   else
 385     dem_mem = (gint16 *)dem_file;
 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 {
 392     g_warning("%s(): file %s does not have right size", __PRETTY_FUNCTION__, basename);
 393     g_mapped_file_unref(mf);
 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++ ) {
 402     dem->n_columns++;
 403     g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
 404     GET_COLUMN(dem,i)->east_west = dem->min_east + arcsec*i;
 405     GET_COLUMN(dem,i)->south = dem->min_north;
 406     GET_COLUMN(dem,i)->n_points = num_rows;
 407     GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*num_rows);
 408   }
 409
 410   int ent = 0;
 411   for ( i = (num_rows - 1); i >= 0; i-- ) {
 412     for ( j = 0; j < num_rows; j++ ) {
 413       GET_COLUMN(dem,j)->points[i] = GINT16_FROM_BE(dem_mem[ent]);
 414       ent++;
 415     }
 416
 417   }
 418
 419   if (zip)
 420     g_free(dem_mem);
 421   g_mapped_file_unref(mf);
 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 {
 429   FILE *f=NULL;
 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
 438   if ( g_access ( file, R_OK ) != 0 )
 439     return NULL;
 440
 441   if ( (strlen(basename)==11 || ((strlen(basename) == 15) && (basename[11] == '.' && basename[12] == 'z' && basename[13] == 'i' && basename[14] == 'p'))) &&
 442        basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
 443        (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W')) {
 444     gboolean is_zip_file = (strlen(basename) == 15);
 445     rv = vik_dem_read_srtm_hgt(file, basename, is_zip_file);
 446     return(rv);
 447   }
 448
 449       /* Create Structure */
 450   rv = g_malloc(sizeof(VikDEM));
 451
 452       /* Header */
 453   f = g_fopen(file, "r");
 454   if ( !f ) {
 455     g_free ( rv );
 456     return NULL;
 457   }
 458   buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
 459   if ( ! dem_parse_header ( buffer, rv ) ) {
 460     g_free ( rv );
 461     fclose(f);
 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);
 490   f = NULL;
 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
 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 */
 554   /* sw */
 555   cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);
 556   rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
 557   ll[0].lon = (dem->min_east + dem->east_scale*cols[0])/3600;
 558   ll[0].lat = (dem->min_north + dem->north_scale*rows[0])/3600;
 559   /* nw */
 560   cols[1] = cols[0];
 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;
 564   /* ne */
 565   cols[2] = cols[0] + 1;
 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;
 569   /* se */
 570   cols[3] = cols[0] + 1;
 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
 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
 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 {
 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_dist*dists[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
 649 }
 650
 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