]>
Commit | Line | Data |
---|---|---|
ad0a8c2d EB |
1 | #include <stdio.h> |
2 | #include <string.h> | |
3 | #include <glib.h> | |
4 | #include <math.h> | |
5 | #include <stdlib.h> | |
46e4fa54 QT |
6 | #include <zlib.h> |
7 | #include <sys/types.h> | |
8 | #include <sys/stat.h> | |
9 | #include <unistd.h> | |
10 | #include <sys/mman.h> | |
11 | ||
ad0a8c2d EB |
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); | |
70147e60 | 90 | /* dem->orig_vert_units = val; now done below */ |
ad0a8c2d EB |
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; | |
70147e60 | 96 | dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS; |
ad0a8c2d EB |
97 | } else { |
98 | dem->east_scale = 3.0; /* arcseconds */ | |
99 | dem->north_scale = 3.0; | |
70147e60 | 100 | dem->orig_vert_units = VIK_DEM_VERT_METERS; |
ad0a8c2d EB |
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 ) { | |
70147e60 EB |
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 | |
ad0a8c2d EB |
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 | ||
46e4fa54 QT |
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) | |
ad0a8c2d | 287 | { |
ad0a8c2d EB |
288 | gint i, j; |
289 | VikDEM *dem; | |
46e4fa54 QT |
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; | |
ad0a8c2d EB |
299 | |
300 | dem = g_malloc(sizeof(VikDEM)); | |
301 | ||
302 | dem->horiz_units = VIK_DEM_HORIZ_LL_ARCSECONDS; | |
70147e60 | 303 | dem->orig_vert_units = VIK_DEM_VERT_DECIMETERS; |
ad0a8c2d EB |
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 | ||
46e4fa54 QT |
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++ ) { | |
ad0a8c2d EB |
357 | dem->n_columns++; |
358 | g_ptr_array_add ( dem->columns, g_malloc(sizeof(VikDEMColumn)) ); | |
46e4fa54 | 359 | GET_COLUMN(dem,i)->east_west = dem->min_east + arcsec*i; |
ad0a8c2d | 360 | GET_COLUMN(dem,i)->south = dem->min_north; |
46e4fa54 QT |
361 | GET_COLUMN(dem,i)->n_points = num_rows; |
362 | GET_COLUMN(dem,i)->points = g_malloc(sizeof(gint16)*num_rows); | |
ad0a8c2d EB |
363 | } |
364 | ||
46e4fa54 QT |
365 | int ent = 0; |
366 | for ( i = (num_rows - 1); i >= 0; i-- ) { | |
367 | for ( j = 0; j < num_rows; j++ ) { | |
792a7488 QT |
368 | GET_COLUMN(dem,j)->points[i] = GINT16_FROM_BE(dem_mem[ent]); |
369 | ent++; | |
ad0a8c2d EB |
370 | } |
371 | ||
372 | } | |
373 | ||
46e4fa54 QT |
374 | if (zip) |
375 | g_free(dem_mem); | |
376 | munmap(dem_file, stat.st_size); | |
ad0a8c2d EB |
377 | return dem; |
378 | } | |
379 | ||
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')) | |
381 | ||
382 | VikDEM *vik_dem_new_from_file(const gchar *file) | |
383 | { | |
384 | FILE *f; | |
385 | VikDEM *rv; | |
386 | gchar buffer[DEM_BLOCK_SIZE+1]; | |
387 | ||
388 | /* use to record state for dem_parse_block */ | |
389 | gint cur_column = -1; | |
390 | gint cur_row = -1; | |
391 | const gchar *basename = a_file_basename(file); | |
392 | ||
393 | /* FILE IO */ | |
394 | f = fopen(file, "r"); | |
395 | if ( !f ) | |
396 | return NULL; | |
397 | ||
46e4fa54 QT |
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' && | |
2dc9ab5e | 400 | (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W')) { |
46e4fa54 QT |
401 | gboolean is_zip_file = (strlen(basename) == 15); |
402 | rv = vik_dem_read_srtm_hgt(f, basename, is_zip_file); | |
2dc9ab5e QT |
403 | fclose(f); |
404 | return(rv); | |
405 | } | |
ad0a8c2d EB |
406 | |
407 | /* Create Structure */ | |
408 | rv = g_malloc(sizeof(VikDEM)); | |
409 | ||
410 | /* Header */ | |
411 | buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0'; | |
412 | if ( ! dem_parse_header ( buffer, rv ) ) { | |
413 | g_free ( rv ); | |
414 | return NULL; | |
415 | } | |
416 | /* TODO: actually use header -- i.e. GET # OF COLUMNS EXPECTED */ | |
417 | ||
418 | rv->columns = g_ptr_array_new(); | |
419 | rv->n_columns = 0; | |
420 | ||
421 | /* Column -- Data */ | |
422 | while (! feof(f) ) { | |
423 | gchar *tmp; | |
424 | ||
425 | /* read block */ | |
426 | buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0'; | |
427 | ||
428 | /* Fix Fortran-style exponentiation */ | |
429 | tmp = buffer; | |
430 | while (*tmp) { | |
431 | if ( *tmp=='D') | |
432 | *tmp='E'; | |
433 | tmp++; | |
434 | } | |
435 | ||
436 | dem_parse_block(buffer, rv, &cur_column, &cur_row); | |
437 | } | |
438 | ||
439 | /* TODO - class C records (right now says 'Invalid' and dies) */ | |
440 | ||
441 | fclose(f); | |
442 | ||
443 | /* 24k scale */ | |
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; | |
446 | ||
447 | /* FIXME bug in 10m DEM's */ | |
448 | if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->north_scale == 10 ) { | |
449 | rv->min_east -= 100; | |
450 | rv->min_north += 200; | |
451 | } | |
452 | ||
453 | ||
454 | return rv; | |
455 | } | |
456 | ||
457 | void vik_dem_free ( VikDEM *dem ) | |
458 | { | |
459 | guint i; | |
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 ); | |
463 | g_free ( dem ); | |
464 | } | |
465 | ||
466 | gint16 vik_dem_get_xy ( VikDEM *dem, guint col, guint row ) | |
467 | { | |
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; | |
472 | } | |
473 | ||
474 | gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north ) | |
475 | { | |
476 | gint col, row; | |
477 | ||
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; | |
481 | ||
482 | col = (gint) floor((east - dem->min_east) / dem->east_scale); | |
483 | row = (gint) floor((north - dem->min_north) / dem->north_scale); | |
484 | ||
485 | return vik_dem_get_xy ( dem, col, row ); | |
486 | } | |
487 | ||
5ef1d57e QT |
488 | static gboolean dem_get_ref_points_elev_dist(VikDEM *dem, |
489 | gdouble east, gdouble north, /* in seconds */ | |
490 | gint16 *elevs, gint16 *dists) | |
491 | { | |
492 | int i; | |
493 | int cols[4], rows[4]; | |
494 | struct LatLon ll[4]; | |
495 | struct LatLon pos; | |
496 | ||
497 | if ( east > dem->max_east || east < dem->min_east || | |
498 | north > dem->max_north || north < dem->min_north ) | |
499 | return FALSE; /* got nothing */ | |
500 | ||
501 | pos.lon = east/3600; | |
502 | pos.lat = north/3600; | |
503 | ||
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 | ||
510 | cols[1] = cols[0]; /*nw*/ | |
511 | rows[1] = rows[0] + 1; | |
512 | ll[1].lon = ll[0].lon; | |
513 | ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600; | |
514 | ||
515 | cols[2] = cols[0] + 1; /*ne*/ | |
516 | rows[2] = rows[0] + 1; | |
517 | ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600; | |
518 | ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600; | |
519 | ||
520 | cols[3] = cols[0] + 1; /*se*/ | |
521 | rows[3] = rows[0]; | |
522 | ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600; | |
523 | ll[3].lat = ll[0].lat; | |
524 | ||
525 | for (i = 0; i < 4; i++) { | |
526 | if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION) | |
527 | return FALSE; | |
528 | dists[i] = a_coords_latlon_diff(&pos, &ll[i]); | |
529 | } | |
530 | ||
531 | return TRUE; /* all OK */ | |
532 | } | |
533 | ||
534 | gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north ) | |
535 | { | |
536 | int i; | |
537 | gint16 elevs[4], dists[4]; | |
538 | ||
539 | if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists)) | |
540 | return VIK_DEM_INVALID_ELEVATION; | |
541 | ||
542 | for (i = 0; i < 4; i++) { | |
543 | if (dists[i] < 1) { | |
544 | return(elevs[i]); | |
545 | } | |
546 | } | |
547 | ||
548 | gdouble t = (gdouble)elevs[0]/dists[0] + (gdouble)elevs[1]/dists[1] + (gdouble)elevs[2]/dists[2] + (gdouble)elevs[3]/dists[3]; | |
549 | gdouble b = 1.0/dists[0] + 1.0/dists[1] + 1.0/dists[2] + 1.0/dists[3]; | |
550 | ||
551 | return(t/b); | |
552 | } | |
553 | ||
554 | gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north ) | |
555 | { | |
556 | return vik_dem_get_simple_interpol(dem, east, north); | |
557 | } | |
558 | ||
ad0a8c2d EB |
559 | void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row ) |
560 | { | |
561 | *col = (gint) floor((east - dem->min_east) / dem->east_scale); | |
562 | *row = (gint) floor((north - dem->min_north) / dem->north_scale); | |
563 | if ( *col < 0 ) *col = 0; | |
564 | if ( *row < 0 ) *row = 0; | |
565 | } | |
566 |