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

   1 /*
   2  * This file is part of mapper
   3  *
   4  * Copyright (C) 2007 Kaj-Michael Lang
   5  *
   6  * Original Algorithms from misc web sites
   7  *
   8  * Default map data provided by http://www.openstreetmap.org/
   9  *
  10  * This program is free software; you can redistribute it and/or modify
  11  * it under the terms of the GNU General Public License as published by
  12  * the Free Software Foundation; either version 2 of the License, or
  13  * (at your option) any later version.
  14  *
  15  * This program is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License along
  21  * with this program; if not, write to the Free Software Foundation, Inc.,
  22  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  23  */
  24
  25 /*
  26  * Different Lat/Lon and related functions
  27  * (conversion, distances, etc)
  28  */
  29 #define _GNU_SOURCE
  30
  31 #include <math.h>
  32 #include <glib.h>
  33
  34 #include "latlon.h"
  35 #include "osm.h"
  36
  37 /* Int ranges for integerized lat/lon */
  38 #define LATLON_MAX 2147483646
  39 #define LATLON_MIN -2147483646
  40
  41 #define EARTH_RADIUS (3440.06479f)
  42
  43 static inline gdouble
  44 magnitude(gdouble x1, gdouble y1, gdouble x2, gdouble y2)
  45 {
  46 gdouble x,y;
  47 x=x2-x1;
  48 y=y2-y1;
  49
  50 return sqrt((x*x)+(y*y));
  51 }
  52
  53 gboolean
  54 distance_point_to_line(gdouble x, gdouble y, gdouble x1, gdouble y1, gdouble x2, gdouble y2, gdouble *d)
  55 {
  56 gdouble lm,u,tmp;
  57 gdouble ix,iy;
  58
  59 lm=magnitude(x1,y1,x2,y2);
  60 if (lm==0.0f)
  61         return FALSE;
  62
  63 tmp=((x-x1)*(x2-x1))+((y-y1)*(y2-y1));
  64 u=tmp/(lm*lm);
  65
  66 if (u<0.0f || u>1.0f)
  67         return FALSE;
  68
  69 ix=x1+u*(x2-x1);
  70 iy=y1+u*(y2-y1);
  71
  72 *d=magnitude(x,y, ix, iy);
  73
  74 return TRUE;
  75 }
  76
  77 /*
  78  * QuadTile calculation functions, from:
  79  * http://trac.openstreetmap.org/browser/sites/rails_port/lib/quad_tile/quad_tile.h
  80  */
  81 guint32
  82 xy2tile(guint x, guint y)
  83 {
  84 guint32 tile=0;
  85 gint i;
  86
  87 for (i=15;i>=0;i--) {
  88         tile=(tile << 1) | ((x >> i) & 1);
  89         tile=(tile << 1) | ((y >> i) & 1);
  90 }
  91 return tile;
  92 }
  93
  94 guint32
  95 lon2x(gdouble lon)
  96 {
  97 return round((lon + 180.0) * 65535.0 / 360.0);
  98 }
  99
 100 guint32
 101 lat2y(gdouble lat)
 102 {
 103 return round((lat + 90.0) * 65535.0 / 180.0);
 104 }
 105
 106 /*
 107  * Convert latitude to integerized+mercator projected value
 108  */
 109 gint32
 110 lat2mp_int(gdouble lat)
 111 {
 112 return lat > 85.051128779 ? LATLON_MAX : lat < -85.051128779 ? LATLON_MIN :
 113         lrint(log(tan(M_PI_4l+lat*M_PIl/360))/M_PIl*LATLON_MAX);
 114 }
 115
 116 /*
 117  * Convert longitude to integerized+mercator projected value
 118  */
 119 gint32
 120 lon2mp_int(gdouble lon)
 121 {
 122 return lrint(lon/180*LATLON_MAX);
 123 }
 124
 125 gdouble
 126 mp_int2lon(gint lon)
 127 {
 128 return (gdouble)lon/LATLON_MAX*180;
 129 }
 130
 131 gdouble
 132 mp_int2lat(gint lat)
 133 {
 134 return 0;
 135 }
 136
 137 /**
 138  * Quick distance calculation, (does not handle earth curvature, so use for quick non exact needs only)
 139  */
 140 gulong
 141 calculate_idistance(gint lat1, gint lon1, gint lat2, gint lon2)
 142 {
 143 return lrint(sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2))));
 144 }
 145
 146 gdouble
 147 calculate_ddistance(gint lat1, gint lon1, gint lat2, gint lon2)
 148 {
 149 return sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2)));
 150 }
 151
 152 /**
 153  * Quick distance for comparing, skips the final square root
 154  */
 155 gulong
 156 calculate_idistance_cmp(gint lat1, gint lon1, gint lat2, gint lon2)
 157 {
 158 return ((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2));
 159 }
 160
 161 /**
 162  * Calculate the distance between two lat/lon pairs.
 163  * The distance is returned in kilometers.
 164  *
 165  */
 166 gdouble
 167 calculate_distance(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
 168 {
 169 gdouble dlat, dlon, slat, slon, a;
 170
 171 /* Convert to radians. */
 172 lat1*=(M_PIl / 180.f);
 173 lon1*=(M_PIl / 180.f);
 174 lat2*=(M_PIl / 180.f);
 175 lon2*=(M_PIl / 180.f);
 176
 177 dlat=lat2 - lat1;
 178 dlon=lon2 - lon1;
 179
 180 slat=sin(dlat / 2.f);
 181 slon=sin(dlon / 2.f);
 182
 183 a=(slat * slat) + (cos(lat1) * cos(lat2) * slon * slon);
 184
 185 return ((2.f * atan2(sqrt(a), sqrt(1.f - a))) * EARTH_RADIUS);
 186 }
 187
 188 /**
 189  * Calculate course from lat1,lon1 to lat2,lon2
 190  *
 191  */
 192 gdouble
 193 calculate_course(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
 194 {
 195 gdouble c,dlon,x,y;
 196
 197 lat1*=(M_PIl / 180.f);
 198 lon1*=(M_PIl / 180.f);
 199 lat2*=(M_PIl / 180.f);
 200 lon2*=(M_PIl / 180.f);
 201
 202 dlon=lon2-lon1;
 203 y=sin(dlon)*cos(lat2);
 204 x=cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon);
 205 c=atan2(y,x);
 206
 207 /* Fixup for compass */
 208 c*=180.f/M_PIl;
 209 c=c+360 % 360;
 210 return c;
 211 }