err.no Git - mapper/blob - src/latlon.c

   1 /*
   2  * Different Lat/Lon related function (conversion, distances, etc)
   3  */
   4 #define _GNU_SOURCE
   5
   6 #include <math.h>
   7 #include <glib.h>
   8
   9 #include "osm.h"
  10
  11 /* Int ranges for integerized lat/lon */
  12 #define LATLON_MAX 2147483646
  13 #define LATLON_MIN -2147483646
  14
  15 #define EARTH_RADIUS (3440.06479f)
  16
  17 /* Convert latitude to integerized+mercator projected value */
  18 gint32
  19 lat2mp_int(gdouble lat)
  20 {
  21 return lat > 85.051128779 ? LATLON_MAX : lat < -85.051128779 ? LATLON_MIN :
  22         lrint(log(tan(M_PI_4l+lat*M_PIl/360))/M_PIl*LATLON_MAX);
  23 }
  24
  25 /* Convert longitude to integerized+mercator projected value */
  26 gint32
  27 lon2mp_int(gdouble lon)
  28 {
  29 return lrint(lon/180*LATLON_MAX);
  30 }
  31
  32 /**
  33  * Calculate the distance between two lat/lon pairs.  The distance is returned
  34  * in kilometer.
  35  */
  36 gfloat calculate_distance(gfloat lat1, gfloat lon1, gfloat lat2, gfloat lon2)
  37 {
  38 gdouble dlat, dlon, slat, slon, a;
  39
  40 /* Convert to radians. */
  41 lat1 *= (M_PI_4l / 180.f);
  42 lon1 *= (M_PI_4l / 180.f);
  43 lat2 *= (M_PI_4l / 180.f);
  44 lon2 *= (M_PI_4l / 180.f);
  45
  46 dlat = lat2 - lat1;
  47 dlon = lon2 - lon1;
  48
  49 slat = sinf(dlat / 2.f);
  50 slon = sinf(dlon / 2.f);
  51 a = (slat * slat) + (cosf(lat1) * cosf(lat2) * slon * slon);
  52
  53 return ((2.f * atan2f(sqrtf(a), sqrtf(1.f - a))) * EARTH_RADIUS);
  54 }
  55