src/dem.c

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