src/dem.c

   1 #include <stdio.h>
   2 #include <string.h>
   3 #include <glib.h>
   4 #include <math.h>
   5 #include <stdlib.h>
   6
   7 #include "dem.h"
   8 #include "file.h"
   9
  10 #define DEM_BLOCK_SIZE 1024
  11 #define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))
  12
  13 static gboolean get_double_and_continue ( gchar **buffer, gdouble *tmp, gboolean warn )
  14 {
  15   gchar *endptr;
  16   *tmp = g_strtod(*buffer, &endptr);
  17   if ( endptr == NULL|| endptr == *buffer ) {
  18     if ( warn )
  19       g_warning("Invalid DEM");
  20     return FALSE;
  21   }
  22   *buffer=endptr;
  23   return TRUE;
  24 }
  25
  26
  27 static gboolean get_int_and_continue ( gchar **buffer, gint *tmp, gboolean warn )
  28 {
  29   gchar *endptr;
  30   *tmp = strtol(*buffer, &endptr, 10);
  31   if ( endptr == NULL|| endptr == *buffer ) {
  32     if ( warn )
  33       g_warning("Invalid DEM");
  34     return FALSE;
  35   }
  36   *buffer=endptr;
  37   return TRUE;
  38 }
  39
  40 static gboolean dem_parse_header ( gchar *buffer, VikDEM *dem )
  41 {
  42   gdouble val;
  43   gint int_val;
  44   guint i;
  45   gchar *tmp = buffer;
  46
  47   /* incomplete header */
  48   if ( strlen(buffer) != 1024 )
  49     return FALSE;
  50
  51   /* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
  52   while (*tmp) {
  53     if ( *tmp=='D')
  54       *tmp='E';
  55     tmp++;
  56   }
  57
  58   /* skip name */
  59   buffer += 149;
  60
  61   /* "DEM level code, pattern code, palaimetric reference system code" -- skip */
  62   get_int_and_continue(&buffer, &int_val, TRUE);
  63   get_int_and_continue(&buffer, &int_val, TRUE);
  64   get_int_and_continue(&buffer, &int_val, TRUE);
  65
  66   /* zone */
  67   get_int_and_continue(&buffer, &int_val, TRUE);
  68   dem->utm_zone = int_val;
  69   /* TODO -- southern or northern hemisphere?! */
  70   dem->utm_letter = 'N';
  71
  72   /* skip numbers 5-19  */
  73   for ( i = 0; i < 15; i++ ) {
  74     if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
  75       g_warning ("Invalid DEM header");
  76       return FALSE;
  77     }
  78   }
  79
  80   /* number 20 -- horizontal unit code (utm/ll) */
  81   get_double_and_continue(&buffer, &val, TRUE);
  82   dem->horiz_units = val;
  83   get_double_and_continue(&buffer, &val, TRUE);
  84   /* dem->orig_vert_units = val; now done below */
  85
  86   /* TODO: do this for real. these are only for 1:24k and 1:250k USGS */
  87   if ( dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS ) {
  88     dem->east_scale = 10.0; /* meters */
  89     dem->north_scale = 10.0;
  90     dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
  91   } else {
  92     dem->east_scale = 3.0; /* arcseconds */
  93     dem->north_scale = 3.0;
  94     dem->orig_vert_units = VIK_DEM_VERT_METERS;
  95   }
  96
  97   /* skip next */
  98   get_double_and_continue(&buffer, &val, TRUE);
  99
 100   /* now we get the four corner points. record the min and max. */
 101   get_double_and_continue(&buffer, &val, TRUE);
 102   dem->min_east = dem->max_east = val;
 103   get_double_and_continue(&buffer, &val, TRUE);
 104   dem->min_north = dem->max_north = val;
 105
 106   for ( i = 0; i < 3; i++ ) {
 107     get_double_and_continue(&buffer, &val, TRUE);
 108     if ( val < dem->min_east ) dem->min_east = val;
 109     if ( val > dem->max_east ) dem->max_east = val;
 110     get_double_and_continue(&buffer, &val, TRUE);
 111     if ( val < dem->min_north ) dem->min_north = val;
 112     if ( val > dem->max_north ) dem->max_north = val;
 113   }
 114
 115   return TRUE;
 116 }
 117
 118 static void dem_parse_block_as_cont ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
 119 {
 120   gint tmp;
 121   while ( *cur_row < GET_COLUMN(dem, *cur_column)->n_points ) {
 122     if ( get_int_and_continue(&buffer, &tmp,FALSE) ) {
 123       if ( dem->orig_vert_units == VIK_DEM_VERT_DECIMETERS )
 124         GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) (tmp / 10);
 125       else
 126         GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) tmp;
 127     } else
 128       return;
 129     (*cur_row)++;
 130   }
 131   *cur_row = -1; /* expecting new column */
 132 }
 133
 134 static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
 135 {
 136   guint n_rows;
 137   gint i;
 138   gdouble east_west, south;
 139   gdouble tmp;
 140
 141   /* 1 x n_rows 1 east_west south x x x DATA */
 142
 143   if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
 144     g_warning("Incorrect DEM Class B record: expected 1");
 145     return;
 146   }
 147
 148   /* don't need this */
 149   if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;
 150
 151   /* n_rows */
 152   if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
 153     return;
 154   n_rows = (guint) tmp;
 155
 156   if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
 157     g_warning("Incorrect DEM Class B record: expected 1");
 158     return;
 159   }
 160
 161   if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
 162     return;
 163   if ( !get_double_and_continue(&buffer, &south, TRUE) )
 164     return;
 165
 166   /* next three things we don't need */
 167   if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
 168   if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
 169   if ( !get_double_and_continue(&buffer, &tmp, TRUE)) return;
 170
 171
 172   dem->n_columns ++;
 173   (*cur_column) ++;
 174
 175   /* empty spaces for things before that were skipped */
 176   (*cur_row) = (south - dem->min_north) / dem->north_scale;
 177   if ( south > dem->max_north || (*cur_row) < 0 )
 178     (*cur_row) = 0;
 179
 180   n_rows += *cur_row;
 181
 182   g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
 183   GET_COLUMN(dem,*cur_column)->east_west = east_west;
 184   GET_COLUMN(dem,*cur_column)->south = south;
 185   GET_COLUMN(dem,*cur_column)->n_points = n_rows;
 186   GET_COLUMN(dem,*cur_column)->points = g_malloc(sizeof(gint16)*n_rows);
 187
 188   /* no information for things before that */
 189   for ( i = 0; i < (*cur_row); i++ )
 190     GET_COLUMN(dem,*cur_column)->points[i] = VIK_DEM_INVALID_ELEVATION;
 191
 192   /* now just continue */
 193   dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );
 194
 195
 196 }
 197
 198 static void dem_parse_block ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
 199 {
 200   /* if haven't read anything or have read all items in a columns and are expecting a new column */
 201   if ( *cur_column == -1 || *cur_row == -1 ) {
 202     dem_parse_block_as_header(buffer, dem, cur_column, cur_row);
 203   } else {
 204     dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
 205   }
 206 }
 207
 208 static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename)
 209 {
 210   gint16 elev;
 211   gint i, j;
 212   VikDEM *dem;
 213
 214
 215   dem = g_malloc(sizeof(VikDEM));
 216
 217   dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
 218   dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
 219   dem->east_scale = dem->north_scale = 3;
 220
 221   /* TODO */
 222   dem->min_north = atoi(basename+1) * 3600;
 223   dem->min_east = atoi(basename+4) * 3600;
 224   if ( basename[0] == 'S' )
 225     dem->min_north = - dem->min_north;
 226   if ( basename[3] == 'W' )
 227     dem->min_east = - dem->min_east;
 228
 229   dem->max_north = 3600 + dem->min_north;
 230   dem->max_east = 3600 + dem->min_east;
 231
 232   dem->columns = g_ptr_array_new();
 233   dem->n_columns = 0;
 234
 235   for ( i = 0; i < 1201; i++ ) {
 236     dem->n_columns++;
 237     g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) );
 238     GET_COLUMN(dem,i)->east_west = dem->min_east + 3*i;
 239     GET_COLUMN(dem,i)->south = dem->min_north;
 240     GET_COLUMN(dem,i)->n_points = 1201;
 241     GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*1201);
 242   }
 243
 244   for ( i = 1200; i >= 0; i-- ) {
 245     for ( j = 0; j < 1201; j++ ) {
 246       if ( feof(f) ) {
 247         g_warning("Incorrect SRTM file: unexpected EOF");
 248         g_print("%d %d\n", i, j);
 249         return dem;
 250       }
 251       fread(&elev, sizeof(elev), 1, f);
 252       gchar *x = (gchar *) &elev;
 253       gchar tmp;
 254       tmp=x[0];
 255       x[0]=x[1];
 256       x[1]=tmp;
 257       GET_COLUMN(dem,j)->points[i] = elev;
 258     }
 259
 260   }
 261
 262   return dem;
 263 }
 264
 265 #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'))
 266
 267 VikDEM *vik_dem_new_from_file(const gchar *file)
 268 {
 269   FILE *f;
 270   VikDEM *rv;
 271   gchar buffer[DEM_BLOCK_SIZE+1];
 272
 273   /* use to record state for dem_parse_block */
 274   gint cur_column = -1;
 275   gint cur_row = -1;
 276   const gchar *basename = a_file_basename(file);
 277
 278       /* FILE IO */
 279   f = fopen(file, "r");
 280   if ( !f )
 281     return NULL;
 282
 283   if ( strlen(basename)==11 && basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
 284        (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W'))
 285     return vik_dem_read_srtm_hgt(f, basename);
 286
 287       /* Create Structure */
 288   rv = g_malloc(sizeof(VikDEM));
 289
 290       /* Header */
 291   buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
 292   if ( ! dem_parse_header ( buffer, rv ) ) {
 293     g_free ( rv );
 294     return NULL;
 295   }
 296   /* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */
 297
 298   rv->columns = g_ptr_array_new();
 299   rv->n_columns = 0;
 300
 301       /* Column -- Data */
 302   while (! feof(f) ) {
 303      gchar *tmp;
 304
 305      /* read block */
 306      buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
 307
 308      /* Fix Fortran-style exponentiation */
 309      tmp = buffer;
 310      while (*tmp) {
 311        if ( *tmp=='D')
 312          *tmp='E';
 313        tmp++;
 314      }
 315
 316      dem_parse_block(buffer, rv, &cur_column, &cur_row);
 317   }
 318
 319      /* TODO - class C records (right now says 'Invalid' and dies) */
 320
 321   fclose(f);
 322
 323   /* 24k scale */
 324   if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
 325     rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;
 326
 327   /* FIXME bug in 10m DEM's */
 328   if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
 329     rv->min_east -= 100;
 330     rv->min_north += 200;
 331   }
 332
 333
 334   return rv;
 335 }
 336
 337 void vik_dem_free ( VikDEM *dem )
 338 {
 339   guint i;
 340   for ( i = 0; i < dem->n_columns; i++)
 341     g_free ( GET_COLUMN(dem, i)->points );
 342   g_ptr_array_free ( dem->columns, TRUE );
 343   g_free ( dem );
 344 }
 345
 346 gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
 347 {
 348   if ( col < dem->n_columns )
 349     if ( row < GET_COLUMN(dem, col)->n_points )
 350       return GET_COLUMN(dem, col)->points[row];
 351   return VIK_DEM_INVALID_ELEVATION;
 352 }
 353
 354 gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
 355 {
 356   gint col, row;
 357
 358   if ( east > dem->max_east || east < dem->min_east ||
 359       north > dem->max_north || north < dem->min_north )
 360     return VIK_DEM_INVALID_ELEVATION;
 361
 362   col = (gint) floor((east - dem->min_east) / dem->east_scale);
 363   row = (gint) floor((north - dem->min_north) / dem->north_scale);
 364
 365   return vik_dem_get_xy ( dem, col, row );
 366 }
 367
 368 void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )
 369 {
 370   *col = (gint) floor((east - dem->min_east) / dem->east_scale);
 371   *row = (gint) floor((north - dem->min_north) / dem->north_scale);
 372   if ( *col < 0 ) *col = 0;
 373   if ( *row < 0 ) *row = 0;
 374 }
 375