]> git.street.me.uk Git - andy/viking.git/blobdiff - src/dem.c
keeping OSM password safe
[andy/viking.git] / src / dem.c
index a3f13ee3e4bd487093bf70968aa1e1f6c02a1f64..9c86a0529b6d476c503ab40498523fa5e2a41a7b 100644 (file)
--- a/src/dem.c
+++ b/src/dem.c
@@ -1,14 +1,48 @@
+/*
+ * viking -- GPS Data and Topo Analyzer, Explorer, and Manager
+ *
+ * Copyright (C) 2003-2008, Evan Battaglia <gtoevan@gmx.net>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <stdio.h>
+#ifdef HAVE_STRING_H
 #include <string.h>
+#endif
 #include <glib.h>
+#ifdef HAVE_MATH_H
 #include <math.h>
+#endif
+#ifdef HAVE_STDLIB_H
 #include <stdlib.h>
+#endif
 #include <zlib.h>
 #include <sys/types.h>
 #include <sys/stat.h>
+#ifdef HAVE_UNISTD_H
 #include <unistd.h>
-#include <sys/mman.h>
+#endif
 
+#include <glib/gstdio.h>
+#include <glib/gi18n.h>
 
 #include "dem.h"
 #include "file.h"
@@ -22,7 +56,7 @@ static gboolean get_double_and_continue ( gchar **buffer, gdouble *tmp, gboolean
   *tmp = g_strtod(*buffer, &endptr);
   if ( endptr == NULL|| endptr == *buffer ) {
     if ( warn )
-      g_warning("Invalid DEM");
+      g_warning(_("Invalid DEM"));
     return FALSE;
   }
   *buffer=endptr;
@@ -36,7 +70,7 @@ static gboolean get_int_and_continue ( gchar **buffer, gint *tmp, gboolean warn
   *tmp = strtol(*buffer, &endptr, 10);
   if ( endptr == NULL|| endptr == *buffer ) {
     if ( warn )
-      g_warning("Invalid DEM");
+      g_warning(_("Invalid DEM"));
     return FALSE;
   }
   *buffer=endptr;
@@ -78,7 +112,7 @@ static gboolean dem_parse_header ( gchar *buffer, VikDEM *dem )
   /* skip numbers 5-19  */
   for ( i = 0; i < 15; i++ ) {
     if ( ! get_double_and_continue(&buffer, &val, FALSE) ) {
-      g_warning ("Invalid DEM header");
+      g_warning (_("Invalid DEM header"));
       return FALSE;
     }
   }
@@ -147,7 +181,7 @@ static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_co
   /* 1 x n_rows 1 east_west south x x x DATA */
 
   if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
-    g_warning("Incorrect DEM Class B record: expected 1");
+    g_warning(_("Incorrect DEM Class B record: expected 1"));
     return;
   }
 
@@ -160,7 +194,7 @@ static void dem_parse_block_as_header ( gchar *buffer, VikDEM *dem, gint *cur_co
   n_rows = (guint) tmp;
 
   if ( (!get_double_and_continue(&buffer, &tmp, TRUE)) || tmp != 1 ) {
-    g_warning("Incorrect DEM Class B record: expected 1");
+    g_warning(_("Incorrect DEM Class B record: expected 1"));
     return;
   }
   
@@ -264,7 +298,7 @@ static void *unzip_hgt_file(gchar *zip_file, gulong *unzip_size)
 
        local_file_header = (struct _lfh *) zip_file;
        if (local_file_header->sig != 0x04034b50) {
-               g_warning("%s(): wrong format\n", __PRETTY_FUNCTION__);
+               g_warning("%s(): wrong format", __PRETTY_FUNCTION__);
                g_free(unzip_data);
                goto end;
        }
@@ -283,19 +317,19 @@ end:
        return(unzip_data);
 }
 
-static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zip)
+static VikDEM *vik_dem_read_srtm_hgt(const gchar *file_name, const gchar *basename, gboolean zip)
 {
   gint i, j;
   VikDEM *dem;
-  struct stat stat;
   off_t file_size;
   gint16 *dem_mem = NULL;
-  gint16 *dem_file = NULL;
+  gchar *dem_file = NULL;
   const gint num_rows_3sec = 1201;
   const gint num_rows_1sec = 3601;
   gint num_rows;
-  int fd = fileno(f);
+  GMappedFile *mf;
   gint arcsec;
+  GError *error = NULL;
 
   dem = g_malloc(sizeof(VikDEM));
 
@@ -316,19 +350,21 @@ static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zi
   dem->columns = g_ptr_array_new();
   dem->n_columns = 0;
 
-  if (fstat(fd, &stat) == -1)
-    g_error("%s(): fstat failed on %s\n", __PRETTY_FUNCTION__, basename);
-  if ((dem_file = mmap(NULL, stat.st_size, PROT_READ, MAP_PRIVATE, fd, 0)) == (void *) -1)
-    g_error("%s(): mmap failed on %s\n", __PRETTY_FUNCTION__, basename);
-
-  file_size = stat.st_size;
-  dem_mem = dem_file;
+  if ((mf = g_mapped_file_new(file_name, FALSE, &error)) == NULL) {
+    g_error(_("Couldn't map file %s: %s"), file_name, error->message);
+    g_error_free(error);
+    g_free(dem);
+    return NULL;
+  }
+  file_size = g_mapped_file_get_length(mf);
+  dem_file = g_mapped_file_get_contents(mf);
+  
   if (zip) {
     void *unzip_mem = NULL;
     gulong ucsize;
 
-    if ((unzip_mem = unzip_hgt_file((gchar *)dem_file, &ucsize)) == NULL) {
-      munmap(dem_file, file_size);
+    if ((unzip_mem = unzip_hgt_file(dem_file, &ucsize)) == NULL) {
+      g_mapped_file_free(mf);
       g_ptr_array_free(dem->columns, TRUE);
       g_free(dem);
       return NULL;
@@ -343,9 +379,8 @@ static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zi
   else if (file_size == (num_rows_1sec * num_rows_1sec * sizeof(gint16)))
     arcsec = 1;
   else {
-    g_warning("%s(): file %s does not have right size\n", __PRETTY_FUNCTION__, basename);
-    munmap(dem_file, file_size);
-    g_ptr_array_free(dem->columns, TRUE);
+    g_warning("%s(): file %s does not have right size", __PRETTY_FUNCTION__, basename);
+    g_mapped_file_free(mf);
     g_free(dem);
     return NULL;
   }
@@ -373,7 +408,7 @@ static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zi
 
   if (zip)
     g_free(dem_mem);
-  munmap(dem_file, stat.st_size);
+  g_mapped_file_free(mf);
   return dem;
 }
 
@@ -381,7 +416,7 @@ static VikDEM *vik_dem_read_srtm_hgt(FILE *f, const gchar *basename, gboolean zi
 
 VikDEM *vik_dem_new_from_file(const gchar *file)
 {
-  FILE *f;
+  FILE *f=NULL;
   VikDEM *rv;
   gchar buffer[DEM_BLOCK_SIZE+1];
 
@@ -390,17 +425,14 @@ VikDEM *vik_dem_new_from_file(const gchar *file)
   gint cur_row = -1;
   const gchar *basename = a_file_basename(file);
 
-      /* FILE IO */
-  f = fopen(file, "r");
-  if ( !f )
+  if ( g_access ( file, R_OK ) != 0 )
     return NULL;
 
   if ( (strlen(basename)==11 || ((strlen(basename) == 15) && (basename[11] == '.' && basename[12] == 'z' && basename[13] == 'i' && basename[14] == 'p'))) &&
        basename[7]=='.' && basename[8]=='h' && basename[9]=='g' && basename[10]=='t' &&
        (basename[0] == 'N' || basename[0] == 'S') && (basename[3] == 'E' || basename[3] =='W')) {
     gboolean is_zip_file = (strlen(basename) == 15);
-    rv = vik_dem_read_srtm_hgt(f, basename, is_zip_file);
-    fclose(f);
+    rv = vik_dem_read_srtm_hgt(file, basename, is_zip_file);
     return(rv);
   }
 
@@ -408,6 +440,11 @@ VikDEM *vik_dem_new_from_file(const gchar *file)
   rv = g_malloc(sizeof(VikDEM));
 
       /* Header */
+  f = g_fopen(file, "r");
+  if ( !f ) {
+    g_free ( rv );
+    return NULL;
+  }
   buffer[fread(buffer, 1, DEM_BLOCK_SIZE, f)] = '\0';
   if ( ! dem_parse_header ( buffer, rv ) ) {
     g_free ( rv );
@@ -439,6 +476,7 @@ VikDEM *vik_dem_new_from_file(const gchar *file)
      /* TODO - class C records (right now says 'Invalid' and dies) */
 
   fclose(f);
+  f = NULL;
 
   /* 24k scale */
   if ( rv->horiz_units == VIK_DEM_HORIZ_UTM_METERS && rv->n_columns >= 2 )
@@ -502,22 +540,23 @@ static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
   pos.lat = north/3600;
 
   /* order of the data: sw, nw, ne, se */
-  cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);  /* sw */
+  /* sw */
+  cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);
   rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
-  ll[0].lon = (east - fabs(fmod(east, dem->east_scale)))/3600;
-  ll[0].lat = (north - fabs(fmod(north, dem->north_scale)))/3600;
-
-  cols[1] = cols[0];         /*nw*/
+  ll[0].lon = (dem->min_east + dem->east_scale*cols[0])/3600;
+  ll[0].lat = (dem->min_north + dem->north_scale*rows[0])/3600;
+  /* nw */
+  cols[1] = cols[0];
   rows[1] = rows[0] + 1;
   ll[1].lon = ll[0].lon;
   ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
-
-  cols[2] = cols[0] + 1;     /*ne*/
+  /* ne */
+  cols[2] = cols[0] + 1;
   rows[2] = rows[0] + 1;
   ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
   ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
-
-  cols[3] = cols[0] + 1;     /*se*/
+  /* se */
+  cols[3] = cols[0] + 1;
   rows[3] = rows[0];
   ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
   ll[3].lat = ll[0].lat;
@@ -528,6 +567,12 @@ static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
     dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
   }
 
+#if 0  /* debug */
+  for (i = 0; i < 4; i++)
+    fprintf(stderr, "%f:%f:%d:%d  ", ll[i].lat, ll[i].lon, dists[i], elevs[i]);
+  fprintf(stderr, "   north_scale=%f\n", dem->north_scale);
+#endif
+
   return TRUE;  /* all OK */
 }
 
@@ -553,7 +598,43 @@ gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
 
 gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
 {
-  return vik_dem_get_simple_interpol(dem, east, north);
+  int i;
+  gint16 elevs[4], dists[4];
+  gint16 max_dist;
+  gdouble t = 0.0;
+  gdouble b = 0.0;
+
+  if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
+    return VIK_DEM_INVALID_ELEVATION;
+
+  max_dist = 0;
+  for (i = 0; i < 4; i++) {
+    if (dists[i] < 1) {
+      return(elevs[i]);
+    }
+    if (dists[i] > max_dist)
+      max_dist = dists[i];
+  }
+
+  gdouble tmp;
+#if 0 /* derived method by Franke & Nielson. Does not seem to work too well here */
+  for (i = 0; i < 4; i++) {
+    tmp = pow((1.0*(max_dist - dists[i])/max_dist*dists[i]), 2);
+    t += tmp*elevs[i];
+    b += tmp;
+  }
+#endif
+
+  for (i = 0; i < 4; i++) {
+    tmp = pow((1.0/dists[i]), 2);
+    t += tmp*elevs[i];
+    b += tmp;
+  }
+
+  // fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);
+
+  return(t/b);
+
 }
 
 void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )