16 #define DEM_BLOCK_SIZE 1024
17 #define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))
19 static gboolean get_double_and_continue ( gchar **buffer, gdouble *tmp, gboolean warn )
22 *tmp = g_strtod(*buffer, &endptr);
23 if ( endptr == NULL|| endptr == *buffer ) {
25 g_warning("Invalid DEM");
33 static gboolean get_int_and_continue ( gchar **buffer, gint *tmp, gboolean warn )
36 *tmp = strtol(*buffer, &endptr, 10);
37 if ( endptr == NULL|| endptr == *buffer ) {
39 g_warning("Invalid DEM");
46 static gboolean dem_parse_header ( gchar *buffer, VikDEM *dem )
53 /* incomplete header */
54 if ( strlen(buffer) != 1024 )
57 /* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
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);
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';
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");
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 */
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;
98 dem->east_scale = 3.0; /* arcseconds */
99 dem->north_scale = 3.0;
100 dem->orig_vert_units = VIK_DEM_VERT_METERS;
104 get_double_and_continue(&buffer, &val, TRUE);
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;
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;
124 static void dem_parse_block_as_cont ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
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);
132 GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) tmp;
137 *cur_row = -1; /* expecting new column */
140 static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
144 gdouble east_west, south;
147 /* 1 x n_rows 1 east_west south x x x DATA */
149 if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
150 g_warning("Incorrect DEM Class B record: expected 1");
154 /* don't need this */
155 if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;
158 if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
160 n_rows = (guint) tmp;
162 if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
163 g_warning("Incorrect DEM Class B record: expected 1");
167 if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
169 if ( !get_double_and_continue(&buffer, &south, TRUE) )
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;
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 )
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);
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;
198 /* now just continue */
199 dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );
204 static void dem_parse_block ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
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);
210 dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
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)
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;
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__);
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);
243 return(stream.total_out);
246 static void *unzip_hgt_file(gchar *zip_file, gulong *unzip_size)
248 void *unzip_data = NULL;
252 guint16 extract_version;
258 guint32 compressed_size;
259 guint32 uncompressed_size;
260 guint16 filename_len;
261 guint16 extra_field_len;
262 } __attribute__ ((__packed__)) *local_file_header = NULL;
265 local_file_header = (struct _lfh *) zip_file;
266 if (local_file_header->sig != 0x04034b50) {
267 g_warning("%s(): wrong format\n", __PRETTY_FUNCTION__);
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;
276 if (!(*unzip_size = uncompress_data(unzip_data, uncompressed_size, zip_data, local_file_header->compressed_size))) {
286 static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zip)
292 gint16 *dem_mem = NULL;
293 gint16 *dem_file = NULL;
294 const gint num_rows_3sec = 1201;
295 const gint num_rows_1sec = 3601;
300 dem = g_malloc(sizeof(VikDEM));
302 dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
303 dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
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;
313 dem->max_north = 3600 + dem->min_north;
314 dem->max_east = 3600 + dem->min_east;
316 dem->columns = g_ptr_array_new();
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);
324 file_size = stat.st_size;
327 void *unzip_mem = NULL;
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);
341 if (file_size == (num_rows_3sec * num_rows_3sec * sizeof(gint16)))
343 else if (file_size == (num_rows_1sec * num_rows_1sec * sizeof(gint16)))
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);
353 num_rows = (arcsec == 3) ? num_rows_3sec : num_rows_1sec;
354 dem->east_scale = dem->north_scale = arcsec;
356 for ( i = 0; i < num_rows; i++ ) {
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);
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]);
376 munmap(dem_file, stat.st_size);
380 #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'))
382 VikDEM *vik_dem_new_from_file(const gchar *file)
386 gchar buffer[DEM_BLOCK_SIZE+1];
388 /* use to record state for dem_parse_block */
389 gint cur_column = -1;
391 const gchar *basename = a_file_basename(file);
394 f = fopen(file, "r");
398 if ( (strlen(basename)==11 || ((strlen(basename) == 15) && (basename[11] == '.' && basename[12] == 'z' && basename[13] == 'i' && basename[14] == 'p'))) &&
399 basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
400 (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W')) {
401 gboolean is_zip_file = (strlen(basename) == 15);
402 rv = vik_dem_read_srtm_hgt(f, basename, is_zip_file);
407 /* Create Structure */
408 rv = g_malloc(sizeof(VikDEM));
411 buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
412 if ( ! dem_parse_header ( buffer, rv ) ) {
416 /* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */
418 rv->columns = g_ptr_array_new();
426 buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
428 /* Fix Fortran-style exponentiation */
436 dem_parse_block(buffer, rv, &cur_column, &cur_row);
439 /* TODO - class C records (right now says 'Invalid' and dies) */
444 if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
445 rv->north_scale = rv->east_scale = GET_COLUMN(rv, 1)->east_west - GET_COLUMN(rv,0)->east_west;
447 /* FIXME bug in 10m DEM's */
448 if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
450 rv->min_north += 200;
457 void vik_dem_free ( VikDEM *dem )
460 for ( i = 0; i < dem->n_columns; i++)
461 g_free ( GET_COLUMN(dem, i)->points );
462 g_ptr_array_free ( dem->columns, TRUE );
466 gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
468 if ( col < dem->n_columns )
469 if ( row < GET_COLUMN(dem, col)->n_points )
470 return GET_COLUMN(dem, col)->points[row];
471 return VIK_DEM_INVALID_ELEVATION;
474 gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
478 if ( east > dem->max_east || east < dem->min_east ||
479 north > dem->max_north || north < dem->min_north )
480 return VIK_DEM_INVALID_ELEVATION;
482 col = (gint) floor((east - dem->min_east) / dem->east_scale);
483 row = (gint) floor((north - dem->min_north) / dem->north_scale);
485 return vik_dem_get_xy ( dem, col, row );
488 static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
489 gdouble east, gdouble north, /* in seconds */
490 gint16 *elevs, gint16 *dists)
493 int cols[4], rows[4];
497 if ( east > dem->max_east || east < dem->min_east ||
498 north > dem->max_north || north < dem->min_north )
499 return FALSE; /* got nothing */
502 pos.lat = north/3600;
504 /* order of the data: sw, nw, ne, se */
505 cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale); /* sw */
506 rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
507 //ll[0].lon = (east - fabs(fmod(east, dem->east_scale)))/3600;
508 //ll[0].lat = (north - fabs(fmod(north, dem->north_scale)))/3600;
509 ll[0].lon = (dem->min_east + dem->east_scale*cols[0])/3600;
510 ll[0].lat = (dem->min_north + dem->north_scale*rows[0])/3600;
512 cols[1] = cols[0]; /*nw*/
513 rows[1] = rows[0] + 1;
514 ll[1].lon = ll[0].lon;
515 ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
517 cols[2] = cols[0] + 1; /*ne*/
518 rows[2] = rows[0] + 1;
519 ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
520 ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
522 cols[3] = cols[0] + 1; /*se*/
524 ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
525 ll[3].lat = ll[0].lat;
527 for (i = 0; i < 4; i++) {
528 if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION)
530 dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
534 for (i = 0; i < 4; i++)
535 fprintf(stderr, "%f:%f:%d:%d ", ll[i].lat, ll[i].lon, dists[i], elevs[i]);
536 fprintf(stderr, " north_scale=%f\n", dem->north_scale);
539 return TRUE; /* all OK */
542 gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
545 gint16 elevs[4], dists[4];
547 if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
548 return VIK_DEM_INVALID_ELEVATION;
550 for (i = 0; i < 4; i++) {
556 gdouble t = (gdouble)elevs[0]/dists[0] + (gdouble)elevs[1]/dists[1] + (gdouble)elevs[2]/dists[2] + (gdouble)elevs[3]/dists[3];
557 gdouble b = 1.0/dists[0] + 1.0/dists[1] + 1.0/dists[2] + 1.0/dists[3];
562 gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
565 gint16 elevs[4], dists[4];
570 if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
571 return VIK_DEM_INVALID_ELEVATION;
574 for (i = 0; i < 4; i++) {
578 if (dists[i] > max_dist)
583 #if 0 /* derived method by Franke & Nielson. Does not seem to work too well here */
584 for (i = 0; i < 4; i++) {
585 tmp = pow((1.0*(max_dist - dists[i])/max_dist*dists[i]), 2);
591 for (i = 0; i < 4; i++) {
592 tmp = pow((1.0/dists[i]), 2);
597 // fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);
603 void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )
605 *col = (gint) floor((east - dem->min_east) / dem->east_scale);
606 *row = (gint) floor((north - dem->min_north) / dem->north_scale);
607 if ( *col < 0 ) *col = 0;
608 if ( *row < 0 ) *row = 0;