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 #include "config.h"
  30
  31 #include <math.h>
  32 #include <glib.h>
  33 #include <glib/gi18n.h>
  34 #include <glib/gstdio.h>
  35
  36 #include "latlon.h"
  37 #include "osm.h"
  38
  39 void
  40 latlon_init(void)
  41 {
  42 DEG_FORMAT_TEXT[DDPDDDDD] = "-dd.ddddd°";
  43 DEG_FORMAT_TEXT[DD_MMPMMM] = "-dd°mm.mmm'";
  44 DEG_FORMAT_TEXT[DD_MM_SSPS] = "-dd°mm'ss.s\"";
  45 DEG_FORMAT_TEXT[DDPDDDDD_NSEW] = "dd.ddddd° S";
  46 DEG_FORMAT_TEXT[DD_MMPMMM_NSEW] = "dd°mm.mmm' S";
  47 DEG_FORMAT_TEXT[DD_MM_SSPS_NSEW] = "dd°mm'ss.s\" S";
  48
  49 UNITS_TEXT[UNITS_KM] = _("km");
  50 UNITS_TEXT[UNITS_MI] = _("mi.");
  51 UNITS_TEXT[UNITS_NM] = _("n.m.");
  52 }
  53
  54 gfloat
  55 angle_diff(gfloat a, gfloat b)
  56 {
  57 gfloat tmax, tmin, t;
  58
  59 tmax=MAX(a,b);
  60 tmin=MIN(a,b);
  61
  62 t=tmax-tmin;
  63 if (t>180) t-=360;
  64
  65 return t;
  66 }
  67
  68 static inline gdouble
  69 magnitude(gdouble x1, gdouble y1, gdouble x2, gdouble y2)
  70 {
  71 gdouble x,y;
  72 x=x2-x1;
  73 y=y2-y1;
  74
  75 return sqrt((x*x)+(y*y));
  76 }
  77
  78 gboolean
  79 distance_point_to_line(gdouble x, gdouble y, gdouble x1, gdouble y1, gdouble x2, gdouble y2, gdouble *d)
  80 {
  81 gdouble lm,u,tmp;
  82 gdouble ix,iy;
  83
  84 lm=magnitude(x1,y1,x2,y2);
  85 if (lm==0.0f)
  86         return FALSE;
  87
  88 tmp=((x-x1)*(x2-x1))+((y-y1)*(y2-y1));
  89 u=tmp/(lm*lm);
  90
  91 if (u<0.0f || u>1.0f)
  92         return FALSE;
  93
  94 ix=x1+u*(x2-x1);
  95 iy=y1+u*(y2-y1);
  96
  97 *d=magnitude(x,y, ix, iy);
  98
  99 return TRUE;
 100 }
 101
 102 /*
 103  * QuadTile calculation functions, from:
 104  * http://trac.openstreetmap.org/browser/sites/rails_port/lib/quad_tile/quad_tile.h
 105  */
 106 guint32
 107 xy2tile(guint x, guint y)
 108 {
 109 guint32 tile=0;
 110 gint i;
 111
 112 for (i=15;i>=0;i--) {
 113         tile=(tile << 1) | ((x >> i) & 1);
 114         tile=(tile << 1) | ((y >> i) & 1);
 115 }
 116 return tile;
 117 }
 118
 119 guint32
 120 lon2x(gdouble lon)
 121 {
 122 return round((lon + 180.0) * 65535.0 / 360.0);
 123 }
 124
 125 guint32
 126 lat2y(gdouble lat)
 127 {
 128 return round((lat + 90.0) * 65535.0 / 180.0);
 129 }
 130
 131 /*
 132  * Convert latitude to integerized+mercator projected value
 133  */
 134 gint32
 135 lat2mp_int(gdouble lat)
 136 {
 137 return lat > 85.051128779 ? LATLON_MAX : lat < -85.051128779 ? LATLON_MIN :
 138         lrint(log(tan(M_PI_4l+lat*M_PIl/360))/M_PIl*LATLON_MAX);
 139 }
 140
 141 /*
 142  * Convert longitude to integerized+mercator projected value
 143  */
 144 gint32
 145 lon2mp_int(gdouble lon)
 146 {
 147 return lrint(lon/180*LATLON_MAX);
 148 }
 149
 150 gdouble
 151 mp_int2lon(gint32 lon)
 152 {
 153 return (gdouble)lon/LATLON_MAX*180;
 154 }
 155
 156 gdouble
 157 mp_int2lat(gint32 lat)
 158 {
 159 return 0;
 160 }
 161
 162 /**
 163  * Quick distance calculation, (does not handle earth curvature, so use for quick non exact needs only)
 164  */
 165 gulong
 166 calculate_idistance(gint lat1, gint lon1, gint lat2, gint lon2)
 167 {
 168 return lrint(sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2))));
 169 }
 170
 171 gdouble
 172 calculate_ddistance(gint lat1, gint lon1, gint lat2, gint lon2)
 173 {
 174 return sqrt((double)((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2)));
 175 }
 176
 177 /**
 178  * Quick distance for comparing, skips the final square root
 179  */
 180 gulong
 181 calculate_idistance_cmp(gint lat1, gint lon1, gint lat2, gint lon2)
 182 {
 183 return ((lat1-lat2)*(lat1-lat2)+(lon1-lon2)*(lon1-lon2));
 184 }
 185
 186 /**
 187  * Calculate the distance between two lat/lon pairs.
 188  * The distance is returned in kilometers.
 189  *
 190  */
 191 gdouble
 192 calculate_distance(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
 193 {
 194 gdouble dlat, dlon, slat, slon, a;
 195
 196 /* Convert to radians. */
 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 dlat=lat2 - lat1;
 203 dlon=lon2 - lon1;
 204
 205 slat=sin(dlat / 2.f);
 206 slon=sin(dlon / 2.f);
 207
 208 a=(slat * slat) + (cos(lat1) * cos(lat2) * slon * slon);
 209
 210 return ((2.f * atan2(sqrt(a), sqrt(1.f - a))) * EARTH_RADIUS);
 211 }
 212
 213 /**
 214  * Calculate course from lat1,lon1 to lat2,lon2
 215  *
 216  */
 217 gdouble
 218 calculate_course_rad(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
 219 {
 220 gdouble c,dlon,x,y;
 221
 222 lat1*=(M_PIl / 180.f);
 223 lon1*=(M_PIl / 180.f);
 224 lat2*=(M_PIl / 180.f);
 225 lon2*=(M_PIl / 180.f);
 226
 227 dlon=lon2-lon1;
 228 y=sin(dlon)*cos(lat2);
 229 x=cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dlon);
 230 c=atan2(y,x);
 231
 232 return c;
 233 }
 234
 235 gdouble
 236 calculate_course(gdouble lat1, gdouble lon1, gdouble lat2, gdouble lon2)
 237 {
 238 gdouble c;
 239
 240 c=calculate_course_rad(lat1, lon1, lat2, lon2);
 241
 242 c*=180.f/M_PIl;
 243 c=c+360 % 360;
 244 return c;
 245 }
 246
 247 void
 248 deg_format(DegFormat degformat, gdouble coor, gchar * scoor, gchar neg_char, gchar pos_char)
 249 {
 250 gdouble min;
 251 gdouble acoor=fabs(coor);
 252
 253 switch (degformat) {
 254         case DDPDDDDD:
 255                 g_sprintf(scoor, "%.5f°", coor);
 256         break;
 257         case DDPDDDDD_NSEW:
 258                 g_sprintf(scoor, "%.5f° %c", acoor, coor < 0.f ? neg_char : pos_char);
 259         break;
 260         case DD_MMPMMM:
 261                 g_sprintf(scoor, "%d°%06.3f'", (int)coor, (acoor - (int)acoor) * 60.0);
 262         break;
 263         case DD_MMPMMM_NSEW:
 264                 g_sprintf(scoor, "%d°%06.3f' %c", (int)acoor, (acoor - (int)acoor) * 60.0,     coor < 0.f ? neg_char : pos_char);
 265         break;
 266         case DD_MM_SSPS:
 267                 min = (acoor - (int)acoor) * 60.0;
 268                 g_sprintf(scoor, "%d°%02d'%04.1f\"", (int)coor, (int)min, ((min - (int)min) * 60.0));
 269         break;
 270         case DD_MM_SSPS_NSEW:
 271                 min = (acoor - (int)acoor) * 60.0;
 272                 g_sprintf(scoor, "%d°%02d'%04.1f\" %c", (int)acoor, (int)min, ((min - (int)min) * 60.0), coor < 0.f ? neg_char : pos_char);
 273         break;
 274         default:
 275                 g_return_if_reached();
 276         break;
 277 }
 278 }
 279