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