X-Git-Url: https://git.street.me.uk/andy/viking.git/blobdiff_plain/5ef1d57eb0a455bf87dc573e6cf3710084630aef..48254848ab0cbe2a9ed0549ecc997c17aa7a925f:/src/dem.c?ds=sidebyside

diff --git a/src/dem.c b/src/dem.c
index a3f13ee3..3f29bce0 100644
--- a/src/dem.c
+++ b/src/dem.c
@@ -10,6 +10,8 @@
 #include <sys/mman.h>
 
 
+#include <glib/gi18n.h>
+
 #include "dem.h"
 #include "file.h"
 
@@ -22,7 +24,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 +38,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 +80,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 +149,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 +162,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 +266,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;
 	}
@@ -343,7 +345,7 @@ 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);
+    g_warning("%s(): file %s does not have right size", __PRETTY_FUNCTION__, basename);
     munmap(dem_file, file_size);
     g_ptr_array_free(dem->columns, TRUE);
     g_free(dem);
@@ -502,22 +504,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 +531,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 +562,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 )