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 gint
  33 calculate_idistance(gint lat1, gint lon1, gint lat2, gint lon2)
  34 {
  35 return lrint(sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2))));
  36 }
  37
  38 /**
  39  * Calculate the distance between two lat/lon pairs.  The distance is returned
  40  * in kilometer.
  41  */
  42 gfloat
  43 calculate_distance(gfloat lat1, gfloat lon1, gfloat lat2, gfloat lon2)
  44 {
  45 gdouble dlat, dlon, slat, slon, a;
  46
  47 /* Convert to radians. */
  48 lat1 *= (M_PI_4l / 180.f);
  49 lon1 *= (M_PI_4l / 180.f);
  50 lat2 *= (M_PI_4l / 180.f);
  51 lon2 *= (M_PI_4l / 180.f);
  52
  53 dlat = lat2 - lat1;
  54 dlon = lon2 - lon1;
  55
  56 slat = sinf(dlat / 2.f);
  57 slon = sinf(dlon / 2.f);
  58 a = (slat * slat) + (cosf(lat1) * cosf(lat2) * slon * slon);
  59
  60 return ((2.f * atan2f(sqrtf(a), sqrtf(1.f - a))) * EARTH_RADIUS);
  61 }
  62