]>
Commit | Line | Data |
---|---|---|
f6171fe3 EB |
1 | /* |
2 | * viking -- GPS Data and Topo Analyzer, Explorer, and Manager | |
3 | * | |
4 | * Copyright (C) 2003-2008, Evan Battaglia <gtoevan@gmx.net> | |
a482007a | 5 | * Copyright (C) 2007, Quy Tonthat <qtonthat@gmail.com> |
f6171fe3 EB |
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 | ||
72b6a3ed JB |
23 | #ifdef HAVE_CONFIG_H |
24 | #include "config.h" | |
25 | #endif | |
26 | ||
ad0a8c2d | 27 | #include <stdio.h> |
8c00358d | 28 | #ifdef HAVE_STRING_H |
ad0a8c2d | 29 | #include <string.h> |
8c00358d | 30 | #endif |
ad0a8c2d | 31 | #include <glib.h> |
8c00358d | 32 | #ifdef HAVE_MATH_H |
ad0a8c2d | 33 | #include <math.h> |
8c00358d GB |
34 | #endif |
35 | #ifdef HAVE_STDLIB_H | |
ad0a8c2d | 36 | #include <stdlib.h> |
8c00358d | 37 | #endif |
46e4fa54 QT |
38 | #include <sys/types.h> |
39 | #include <sys/stat.h> | |
31766324 MA |
40 | #ifdef HAVE_UNISTD_H |
41 | #include <unistd.h> | |
42 | #endif | |
46e4fa54 | 43 | |
6af463da | 44 | #include <glib/gstdio.h> |
4c77d5e0 GB |
45 | #include <glib/gi18n.h> |
46 | ||
ba048e02 | 47 | #include "compression.h" |
ad0a8c2d | 48 | #include "dem.h" |
a660c6aa RN |
49 | #include "coords.h" |
50 | #include "fileutils.h" | |
ad0a8c2d | 51 | |
316945d8 GB |
52 | /* Compatibility */ |
53 | #if ! GLIB_CHECK_VERSION(2,22,0) | |
54 | #define g_mapped_file_unref g_mapped_file_free | |
55 | #endif | |
56 | ||
ad0a8c2d EB |
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 ) | |
4c77d5e0 | 66 | g_warning(_("Invalid DEM")); |
ad0a8c2d EB |
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 ) | |
4c77d5e0 | 80 | g_warning(_("Invalid DEM")); |
ad0a8c2d EB |
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) ) { | |
4c77d5e0 | 122 | g_warning (_("Invalid DEM header")); |
ad0a8c2d EB |
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); | |
70147e60 | 131 | /* dem->orig_vert_units = val; now done below */ |
ad0a8c2d EB |
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; | |
70147e60 | 137 | dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS; |
ad0a8c2d EB |
138 | } else { |
139 | dem->east_scale = 3.0; /* arcseconds */ | |
140 | dem->north_scale = 3.0; | |
70147e60 | 141 | dem->orig_vert_units = VIK_DEM_VERT_METERS; |
ad0a8c2d EB |
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 ) { | |
70147e60 EB |
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 | |
ad0a8c2d EB |
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 ) { | |
4c77d5e0 | 191 | g_warning(_("Incorrect DEM Class B record: expected 1")); |
ad0a8c2d EB |
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 ) { | |
4c77d5e0 | 204 | g_warning(_("Incorrect DEM Class B record: expected 1")); |
ad0a8c2d EB |
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 | ||
31766324 | 255 | static VikDEM *vik_dem_read_srtm_hgt(const gchar *file_name, const gchar *basename, gboolean zip) |
ad0a8c2d | 256 | { |
ad0a8c2d EB |
257 | gint i, j; |
258 | VikDEM *dem; | |
46e4fa54 QT |
259 | off_t file_size; |
260 | gint16 *dem_mem = NULL; | |
31766324 | 261 | gchar *dem_file = NULL; |
46e4fa54 QT |
262 | const gint num_rows_3sec = 1201; |
263 | const gint num_rows_1sec = 3601; | |
264 | gint num_rows; | |
31766324 | 265 | GMappedFile *mf; |
46e4fa54 | 266 | gint arcsec; |
a95be497 | 267 | GError *error = NULL; |
ad0a8c2d EB |
268 | |
269 | dem = g_malloc(sizeof(VikDEM)); | |
270 | ||
271 | dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS; | |
70147e60 | 272 | dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS; |
ad0a8c2d EB |
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 | ||
a95be497 | 288 | if ((mf = g_mapped_file_new(file_name, FALSE, &error)) == NULL) { |
24ff9c7b | 289 | g_critical(_("Couldn't map file %s: %s"), file_name, error->message); |
a95be497 | 290 | g_error_free(error); |
31766324 MA |
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 | ||
46e4fa54 QT |
297 | if (zip) { |
298 | void *unzip_mem = NULL; | |
299 | gulong ucsize; | |
300 | ||
ba048e02 | 301 | if ((unzip_mem = unzip_file(dem_file, &ucsize)) == NULL) { |
316945d8 | 302 | g_mapped_file_unref(mf); |
9ddb3bd9 | 303 | g_ptr_array_foreach ( dem->columns, (GFunc)g_free, NULL ); |
46e4fa54 QT |
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 | } | |
92390c64 CK |
312 | else |
313 | dem_mem = (gint16 *)dem_file; | |
46e4fa54 QT |
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 { | |
4258f4e2 | 320 | g_warning("%s(): file %s does not have right size", __PRETTY_FUNCTION__, basename); |
316945d8 | 321 | g_mapped_file_unref(mf); |
46e4fa54 QT |
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++ ) { | |
ad0a8c2d EB |
330 | dem->n_columns++; |
331 | g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) ); | |
46e4fa54 | 332 | GET_COLUMN(dem,i)->east_west = dem->min_east + arcsec*i; |
ad0a8c2d | 333 | GET_COLUMN(dem,i)->south = dem->min_north; |
46e4fa54 QT |
334 | GET_COLUMN(dem,i)->n_points = num_rows; |
335 | GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*num_rows); | |
ad0a8c2d EB |
336 | } |
337 | ||
46e4fa54 QT |
338 | int ent = 0; |
339 | for ( i = (num_rows - 1); i >= 0; i-- ) { | |
340 | for ( j = 0; j < num_rows; j++ ) { | |
792a7488 QT |
341 | GET_COLUMN(dem,j)->points[i] = GINT16_FROM_BE(dem_mem[ent]); |
342 | ent++; | |
ad0a8c2d EB |
343 | } |
344 | ||
345 | } | |
346 | ||
46e4fa54 QT |
347 | if (zip) |
348 | g_free(dem_mem); | |
316945d8 | 349 | g_mapped_file_unref(mf); |
ad0a8c2d EB |
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 | { | |
31766324 | 357 | FILE *f=NULL; |
ad0a8c2d EB |
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 | ||
31766324 | 366 | if ( g_access ( file, R_OK ) != 0 ) |
ad0a8c2d EB |
367 | return NULL; |
368 | ||
46e4fa54 QT |
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' && | |
2dc9ab5e | 371 | (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W')) { |
46e4fa54 | 372 | gboolean is_zip_file = (strlen(basename) == 15); |
31766324 | 373 | rv = vik_dem_read_srtm_hgt(file, basename, is_zip_file); |
2dc9ab5e QT |
374 | return(rv); |
375 | } | |
ad0a8c2d EB |
376 | |
377 | /* Create Structure */ | |
378 | rv = g_malloc(sizeof(VikDEM)); | |
379 | ||
380 | /* Header */ | |
31766324 | 381 | f = g_fopen(file, "r"); |
3b36279c JJ |
382 | if ( !f ) { |
383 | g_free ( rv ); | |
31766324 | 384 | return NULL; |
3b36279c | 385 | } |
ad0a8c2d EB |
386 | buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0'; |
387 | if ( ! dem_parse_header ( buffer, rv ) ) { | |
388 | g_free ( rv ); | |
dad45413 | 389 | fclose(f); |
ad0a8c2d EB |
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); | |
8c060406 | 418 | f = NULL; |
ad0a8c2d EB |
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 ); | |
9ddb3bd9 | 439 | g_ptr_array_foreach ( dem->columns, (GFunc)g_free, NULL ); |
ad0a8c2d EB |
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 | ||
5ef1d57e QT |
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 */ | |
245575e1 QT |
483 | /* sw */ |
484 | cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale); | |
5ef1d57e | 485 | rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale); |
228213c5 QT |
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; | |
245575e1 QT |
488 | /* nw */ |
489 | cols[1] = cols[0]; | |
5ef1d57e QT |
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; | |
245575e1 QT |
493 | /* ne */ |
494 | cols[2] = cols[0] + 1; | |
5ef1d57e QT |
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; | |
245575e1 QT |
498 | /* se */ |
499 | cols[3] = cols[0] + 1; | |
5ef1d57e QT |
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 | ||
228213c5 QT |
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 | ||
5ef1d57e QT |
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 | { | |
228213c5 QT |
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 | ||
5ef1d57e QT |
578 | } |
579 | ||
ad0a8c2d EB |
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 |