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 gdouble
  33 mp_int2lon(gint lon)
  34 {
  35 return (gdouble)lon/LATLON_MAX*180;
  36 }
  37
  38 gdouble
  39 mp_int2lat(gint lat)
  40 {
  41
  42 }
  43
  44 gint
  45 calculate_idistance(gint lat1, gint lon1, gint lat2, gint lon2)
  46 {
  47 return lrint(sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2))));
  48 }
  49
  50 /**
  51  * Quick distance for comparing, skips the final square root
  52  */
  53 gint
  54 calculate_idistance_cmp(gint lat1, gint lon1, gint lat2, gint lon2)
  55 {
  56 return ((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2));
  57 }
  58
  59 /**
  60  * Calculate the distance between two lat/lon pairs.
  61  * The distance is returned in kilometers.
  62  *
  63  */
  64 gdouble
  65 calculate_distance(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
  66 {
  67 gdouble dlat, dlon, slat, slon, a;
  68
  69 /* Convert to radians. */
  70 lat1*=(M_PI_4l / 180.f);
  71 lon1*=(M_PI_4l / 180.f);
  72 lat2*=(M_PI_4l / 180.f);
  73 lon2*=(M_PI_4l / 180.f);
  74
  75 dlat=lat2 - lat1;
  76 dlon=lon2 - lon1;
  77
  78 slat=sinf(dlat / 2.f);
  79 slon=sinf(dlon / 2.f);
  80
  81 a=(slat * slat) + (cosf(lat1) * cosf(lat2) * slon * slon);
  82
  83 return ((2.f * atan2f(sqrtf(a), sqrtf(1.f - a))) * EARTH_RADIUS);
  84 }
  85