2 * viking -- GPS Data and Topo Analyzer, Explorer, and Manager
4 * Copyright (C) 2003-2008, Evan Battaglia <gtoevan@gmx.net>
5 * Copyright (C) 2007, Quy Tonthat <qtonthat@gmail.com>
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.
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.
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
38 #include <sys/types.h>
44 #include <glib/gstdio.h>
45 #include <glib/gi18n.h>
47 #include "compression.h"
50 #include "fileutils.h"
53 #if ! GLIB_CHECK_VERSION(2,22,0)
54 #define g_mapped_file_unref g_mapped_file_free
57 #define DEM_BLOCK_SIZE 1024
58 #define GET_COLUMN(dem,n) ((VikDEMColumn *)g_ptr_array_index( (dem)->columns, (n) ))
60 static gboolean get_double_and_continue ( gchar **buffer, gdouble *tmp, gboolean warn )
63 *tmp = g_strtod(*buffer, &endptr);
64 if ( endptr == NULL|| endptr == *buffer ) {
66 g_warning(_("Invalid DEM"));
74 static gboolean get_int_and_continue ( gchar **buffer, gint *tmp, gboolean warn )
77 *tmp = strtol(*buffer, &endptr, 10);
78 if ( endptr == NULL|| endptr == *buffer ) {
80 g_warning(_("Invalid DEM"));
87 static gboolean dem_parse_header ( gchar *buffer, VikDEM *dem )
94 /* incomplete header */
95 if ( strlen(buffer) != 1024 )
98 /* fix Fortran-style exponentiation 1.0D5 -> 1.0E5 */
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);
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';
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"));
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 */
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;
139 dem->east_scale = 3.0; /* arcseconds */
140 dem->north_scale = 3.0;
141 dem->orig_vert_units = VIK_DEM_VERT_METERS;
145 get_double_and_continue(&buffer, &val, TRUE);
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;
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;
165 static void dem_parse_block_as_cont ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
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);
173 GET_COLUMN(dem, *cur_column)->points[*cur_row] = (gint16) tmp;
178 *cur_row = -1; /* expecting new column */
181 static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
185 gdouble east_west, south;
188 /* 1 x n_rows 1 east_west south x x x DATA */
190 if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
191 g_warning(_("Incorrect DEM Class B record: expected 1"));
195 /* don't need this */
196 if ( !get_double_and_continue(&buffer, &tmp, TRUE ) ) return;
199 if ( !get_double_and_continue(&buffer, &tmp, TRUE ) )
201 n_rows = (guint) tmp;
203 if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
204 g_warning(_("Incorrect DEM Class B record: expected 1"));
208 if ( !get_double_and_continue(&buffer, &east_west, TRUE) )
210 if ( !get_double_and_continue(&buffer, &south, TRUE) )
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;
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 )
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);
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;
239 /* now just continue */
240 dem_parse_block_as_cont ( buffer, dem, cur_column, cur_row );
245 static void dem_parse_block ( gchar *buffer, VikDEM *dem, gint *cur_column, gint *cur_row )
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);
251 dem_parse_block_as_cont(buffer, dem, cur_column, cur_row);
255 static VikDEM *vik_dem_read_srtm_hgt(const gchar *file_name, const gchar *basename, gboolean zip)
260 gint16 *dem_mem = NULL;
261 gchar *dem_file = NULL;
262 const gint num_rows_3sec = 1201;
263 const gint num_rows_1sec = 3601;
267 GError *error = NULL;
269 dem = g_malloc(sizeof(VikDEM));
271 dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS;
272 dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS;
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;
282 dem->max_north = 3600 + dem->min_north;
283 dem->max_east = 3600 + dem->min_east;
285 dem->columns = g_ptr_array_new();
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);
294 file_size = g_mapped_file_get_length(mf);
295 dem_file = g_mapped_file_get_contents(mf);
298 void *unzip_mem = NULL;
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);
313 dem_mem = (gint16 *)dem_file;
315 if (file_size == (num_rows_3sec * num_rows_3sec * sizeof(gint16)))
317 else if (file_size == (num_rows_1sec * num_rows_1sec * sizeof(gint16)))
320 g_warning("%s(): file %s does not have right size", __PRETTY_FUNCTION__, basename);
321 g_mapped_file_unref(mf);
326 num_rows = (arcsec == 3) ? num_rows_3sec : num_rows_1sec;
327 dem->east_scale = dem->north_scale = arcsec;
329 for ( i = 0; i < num_rows; i++ ) {
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);
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]);
349 g_mapped_file_unref(mf);
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'))
355 VikDEM *vik_dem_new_from_file(const gchar *file)
359 gchar buffer[DEM_BLOCK_SIZE+1];
361 /* use to record state for dem_parse_block */
362 gint cur_column = -1;
364 const gchar *basename = a_file_basename(file);
366 if ( g_access ( file, R_OK ) != 0 )
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);
377 /* Create Structure */
378 rv = g_malloc(sizeof(VikDEM));
381 f = g_fopen(file, "r");
386 buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
387 if ( ! dem_parse_header ( buffer, rv ) ) {
392 /* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */
394 rv->columns = g_ptr_array_new();
402 buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
404 /* Fix Fortran-style exponentiation */
412 dem_parse_block(buffer, rv, &cur_column, &cur_row);
415 /* TODO - class C records (right now says 'Invalid' and dies) */
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;
424 /* FIXME bug in 10m DEM's */
425 if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) {
427 rv->min_north += 200;
434 void vik_dem_free ( VikDEM *dem )
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 );
444 gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row )
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;
452 gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
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;
460 col = (gint) floor((east - dem->min_east) / dem->east_scale);
461 row = (gint) floor((north - dem->min_north) / dem->north_scale);
463 return vik_dem_get_xy ( dem, col, row );
466 static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
467 gdouble east, gdouble north, /* in seconds */
468 gint16 *elevs, gint16 *dists)
471 int cols[4], rows[4];
475 if ( east > dem->max_east || east < dem->min_east ||
476 north > dem->max_north || north < dem->min_north )
477 return FALSE; /* got nothing */
480 pos.lat = north/3600;
482 /* order of the data: sw, nw, ne, se */
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;
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;
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;
499 cols[3] = cols[0] + 1;
501 ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
502 ll[3].lat = ll[0].lat;
504 for (i = 0; i < 4; i++) {
505 if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION)
507 dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
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);
516 return TRUE; /* all OK */
519 gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
522 gint16 elevs[4], dists[4];
524 if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
525 return VIK_DEM_INVALID_ELEVATION;
527 for (i = 0; i < 4; i++) {
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];
539 gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
542 gint16 elevs[4], dists[4];
547 if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
548 return VIK_DEM_INVALID_ELEVATION;
551 for (i = 0; i < 4; i++) {
555 if (dists[i] > max_dist)
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);
568 for (i = 0; i < 4; i++) {
569 tmp = pow((1.0/dists[i]), 2);
574 // fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);
580 void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )
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;