Pileus Git - grits/blob - src/grits-util.c

   1 /*
   2  * Copyright (C) 2009-2011 Andy Spencer <andy753421@gmail.com>
   3  *
   4  * This program is free software: you can redistribute it and/or modify
   5  * it under the terms of the GNU General Public License as published by
   6  * the Free Software Foundation, either version 3 of the License, or
   7  * (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  16  */
  17
  18 /**
  19  * SECTION:grits-util
  20  * @short_description: Geographic utilities
  21  *
  22  * Miscellaneous utility functions, these deal mostly with coordinate
  23  * conversion. Below are some examples that should help demonstrate how these
  24  * functions work.
  25  *
  26  * <example>
  27  * <title>Terminology</title>
  28  * <programlisting>
  29  * deg    - Degrees
  30  * rad    - Radians, also radius
  31  * m      - Meters, for earth-based distances
  32  * px     - Pixels, for screen-based distances
  33  *
  34  * height - Height, the distance above the geoid (ground)
  35  * elev   - Elevation, the distance above the spheroid
  36  * rad    - Radius, the distance from the center of the earth
  37  *
  38  * lat    - Latitude, amount north-south, -90 (S) .. 90 (N)
  39  * lon    - Longitude, amount east-west, -180 (W) .. 180 (E)
  40  * incl   - Inclination, polar equiv of  latitude, Pi .. 0
  41  * azim   - Azimuth, polar equiv of longitude, -Pi .. Pi
  42  *
  43  * x      -  0° lon is positive
  44  * y      - 90° lon is positive
  45  * z      - North pole is positive
  46  *
  47  * llh    - lat,lon,height
  48  * lle    - lat,lon,elev
  49  * llr    - lat,lon,rad
  50  * pol    - incl,azim,rad
  51  * xyz    - x,y,z
  52  * </programlisting>
  53  * </example>
  54  *
  55  * <example>
  56  * <title>Conversions</title>
  57  * <programlisting>
  58  *             lat    lon   elev ->      x      y      z
  59  * lle2xyz:    0.0,   0.0,   0.0 ->    0.0,   0.0,  10.0
  60  * lle2xyz:   90.0,   0.0,   0.0 ->    0.0,  10.0,   0.0
  61  * lle2xyz:    0.0,  90.0,   0.0 ->   10.0,   0.0,   0.0
  62  *
  63  *               x      y      z ->    lat    lon   elev
  64  * xyz2lle:   10.0,   0.0,   0.0 ->    0.0,  90.0,   0.0
  65  * xyz2lle:    0.0,  10.0,   0.0 ->   90.0,   0.0,   0.0
  66  * xyz2lle:    0.0,   0.0,  10.0 ->    0.0,   0.0,   0.0
  67  * </programlisting>
  68  * </example>
  69  */
  70
  71 #include <glib.h>
  72 #include <math.h>
  73 #include <stdio.h>
  74
  75 #include "grits-util.h"
  76
  77 /**************
  78  * GritsPoint *
  79  **************/
  80 /**
  81  * grits_point_set_lle:
  82  * @point: the point to modify
  83  * @lat:   the new latitude
  84  * @lon:   the new longitude
  85  * @elev:  the new elevation
  86  *
  87  * Set the latitude, longitude, and elevation for a point.
  88  */
  89 void grits_point_set_lle(GritsPoint *point, gdouble lat, gdouble lon, gdouble elev)
  90 {
  91         point->lat  = lat;
  92         point->lon  = lon;
  93         point->elev = elev;
  94 }
  95
  96
  97 /***************
  98  * GritsBounds *
  99  ***************/
 100 /**
 101  * grits_bounds_set_bounds:
 102  * @n: the north edge
 103  * @s: the south edge
 104  * @e: the east edge
 105  * @w: the west edge
 106  *
 107  * Set the north, south, east, and west edges of the bounding box
 108  */
 109 void grits_bounds_set_bounds(GritsBounds *bounds,
 110                 gdouble n, gdouble s, gdouble e, gdouble w)
 111 {
 112         bounds->n = n;
 113         bounds->s = s;
 114         bounds->e = e;
 115         bounds->w = w;
 116 }
 117
 118
 119 /******************
 120  * Global helpers *
 121  ******************/
 122 /**
 123  * lle2xyz:
 124  * @lat:  the latitude
 125  * @lon:  the longitude
 126  * @elev: the elevation
 127  * @x:    the resulting x coordinate
 128  * @y:    the resulting y coordinate
 129  * @z:    the resulting z coordinate
 130  *
 131  * Convert a point from latitude, longitude, and elevation to x, y and z
 132  * coordinates.
 133  */
 134 void lle2xyz(gdouble lat, gdouble lon, gdouble elev,
 135                 gdouble *x, gdouble *y, gdouble *z)
 136 {
 137         gdouble rad  = elev2rad(elev);
 138         gdouble azim = lon2azim(lon);
 139         gdouble incl = lat2incl(lat);
 140         *z = rad * cos(azim) * sin(incl);
 141         *x = rad * sin(azim) * sin(incl);
 142         *y = rad * cos(incl);
 143 }
 144
 145 /**
 146  * xyz2lle:
 147  * @x:    the x coordinate
 148  * @y:    the y coordinate
 149  * @z:    the z coordinate
 150  * @lat:  the resulting latitude
 151  * @lon:  the resulting longitude
 152  * @elev: the resulting elevation
 153  *
 154  * Convert a point from x, y and z coordinates to latitude, longitude, and
 155  * elevation.
 156  */
 157 void xyz2lle(gdouble x, gdouble y, gdouble z,
 158                 gdouble *lat, gdouble *lon, gdouble *elev)
 159 {
 160         gdouble rad = sqrt(x*x + y*y + z*z);
 161         *lat  = incl2lat(acos(y / rad));
 162         *lon  = azim2lon(atan2(x,z));
 163         *elev = rad2elev(rad);
 164 }
 165
 166 /**
 167  * xyz2ll:
 168  * @x:    the x coordinate
 169  * @y:    the y coordinate
 170  * @z:    the z coordinate
 171  * @lat:  the resulting latitude
 172  * @lon:  the resulting longitude
 173  *
 174  * Get the latitude and longitude for a x, y, z value.
 175  */
 176 void xyz2ll(gdouble x, gdouble y, gdouble z,
 177                 gdouble *lat, gdouble *lon)
 178 {
 179         gdouble rad = sqrt(x*x + y*y + z*z);
 180         *lat = incl2lat(acos(y / rad));
 181         *lon = azim2lon(atan2(x,z));
 182 }
 183
 184 /**
 185  * ll2m:
 186  * @lon_dist: the distance in degrees of longitude
 187  * @lat:      the latitude to calculate at
 188  *
 189  * Calculate the distance of longitudinal span at a particular latitude.
 190  *
 191  * Returns: the distance in meters
 192  */
 193 gdouble ll2m(gdouble lon_dist, gdouble lat)
 194 {
 195         gdouble incl = lat2incl(lat);
 196         gdouble rad  = sin(incl) * EARTH_R;
 197         gdouble circ = 2 * G_PI * rad;
 198         return lon_dist/360 * circ;
 199 }
 200
 201 /**
 202  * distd:
 203  * @a: the first point
 204  * @b: the second point
 205  *
 206  * Calculate the distance between two three dimensional points.
 207  *
 208  * Returns: the distance between the points
 209  */
 210 gdouble distd(gdouble *a, gdouble *b)
 211 {
 212         return sqrt((a[0]-b[0])*(a[0]-b[0]) +
 213                     (a[1]-b[1])*(a[1]-b[1]) +
 214                     (a[2]-b[2])*(a[2]-b[2]));
 215 }
 216
 217 /**
 218  * lon_avg:
 219  * @a: the first longitude
 220  * @b: the second longitude
 221  *
 222  * Calculate the average longitude between two longitudes. This is smart about
 223  * which side of the globe the resulting longitude is placed on.
 224  *
 225  * Returns: the average
 226  */
 227 gdouble lon_avg(gdouble a, gdouble b)
 228 {
 229         gdouble diff = ABS(a-b);
 230         gdouble avg  = (a+b)/2;
 231         if (diff > 180) {
 232                 if (avg >= 0)
 233                         avg -= 180;
 234                 else
 235                         avg += 180;
 236         }
 237         return avg;
 238 }
 239
 240 /**
 241  * crossd3:
 242  * @a:   the left   point
 243  * @b:   the center point
 244  * @c:   the right  point
 245  * @out: the cross product
 246  *
 247  * Calculate the cross product of three points
 248  */
 249 void crossd3(gdouble *a, gdouble *b, gdouble *c, gdouble *out)
 250 {
 251         double ba[3], bc[3];
 252
 253         ba[0] = a[0] - b[0];
 254         ba[1] = a[1] - b[1];
 255         ba[2] = a[2] - b[2];
 256
 257         bc[0] = c[0] - b[0];
 258         bc[1] = c[1] - b[1];
 259         bc[2] = c[2] - b[2];
 260
 261         crossd(ba, bc, out);
 262 }
 263
 264 /**
 265  * crossd:
 266  * @a: the first vector
 267  * @b: the second vector
 268  * @c: the cross product
 269  *
 270  * Calculate the cross product of two vectors
 271  */
 272 void crossd(gdouble *a, gdouble *b, gdouble *out)
 273 {
 274         out[0] = a[1] * b[2] - a[2] * b[1];
 275         out[1] = a[2] * b[0] - a[0] * b[2];
 276         out[2] = a[0] * b[1] - a[1] * b[0];
 277 }
 278
 279 /**
 280  * lengthd:
 281  * @a: the vector
 282  *
 283  * Calculate the length (magnitude) of a vector.
 284  *
 285  * Returns: the length
 286  */
 287 gdouble lengthd(gdouble *a)
 288 {
 289         return sqrt(a[0] * a[0] +
 290                     a[1] * a[1] +
 291                     a[2] * a[2]);
 292 }
 293
 294 /**
 295  * normd:
 296  * @a: the first longitude
 297  *
 298  * Normalize a vector to a mangitude of 1
 299  */
 300 void normd(gdouble *a)
 301 {
 302         gdouble total = lengthd(a);
 303         a[0] = a[0] / total;
 304         a[1] = a[1] / total;
 305         a[2] = a[2] / total;
 306 }
 307
 308 /**
 309  * parse_points:
 310  * @string:    String representation of the points
 311  * @group_sep: Group separator
 312  * @point_sep: Point separator
 313  * @coord_sep: Coordinate separator
 314  * @bounds:    The bounding box of all the points, or NULL
 315  * @center:    The center of the @bounds, or NULL
 316  *
 317  * Parse a string of the form:
 318  *   string -> group [@group_sep group] ...
 319  *   group  -> point [@point_sep point] ...
 320  *   point  -> latitude @coord_sep longitude [@coord_sep elevation]
 321  *
 322  * For example
 323  *   parse_points("30,-80 30,-120 50,-120 50,-80", "\t", " ", ",");
 324  *
 325  * Returns: zero-terminated array of groups of points
 326  */
 327 GritsPoints *parse_points(const gchar *string,
 328                 const gchar *group_sep, const gchar *point_sep, const gchar *coord_sep,
 329                 GritsBounds *bounds, GritsPoint *center)
 330 {
 331         /* Split and count groups */
 332         gchar **sgroups = g_strsplit(string, group_sep, -1);
 333         int     ngroups = g_strv_length(sgroups);
 334
 335         GritsBounds _bounds = {-90, 90, -180, 180};
 336         gdouble (**groups)[3] = (gpointer)g_new0(double*, ngroups+1);
 337         for (int pi = 0; pi < ngroups; pi++) {
 338                 /* Split and count coordinates */
 339                 gchar **scoords = g_strsplit(sgroups[pi], point_sep, -1);
 340                 int     ncoords = g_strv_length(scoords);
 341
 342                 /* Create binary coords */
 343                 gdouble (*coords)[3] = (gpointer)g_new0(gdouble, 3*(ncoords+1));
 344                 for (int ci = 0; ci < ncoords; ci++) {
 345                         gdouble lat, lon;
 346                         sscanf(scoords[ci], "%lf,%lf", &lat, &lon);
 347                         if (bounds || center) {
 348                                 if (lat > _bounds.n) _bounds.n = lat;
 349                                 if (lat < _bounds.s) _bounds.s = lat;
 350                                 if (lon > _bounds.e) _bounds.e = lon;
 351                                 if (lon < _bounds.w) _bounds.w = lon;
 352                         }
 353                         lle2xyz(lat, lon, 0,
 354                                         &coords[ci][0],
 355                                         &coords[ci][1],
 356                                         &coords[ci][2]);
 357                 }
 358
 359                 /* Insert coords into line array */
 360                 groups[pi] = coords;
 361                 g_strfreev(scoords);
 362         }
 363         g_strfreev(sgroups);
 364
 365         /* Output */
 366         if (bounds)
 367                 *bounds = _bounds;
 368         if (center) {
 369                 center->lat  = (_bounds.n + _bounds.s)/2;
 370                 center->lon  = lon_avg(_bounds.e, _bounds.w);
 371                 center->elev = 0;
 372         }
 373         return groups;
 374 }
 375
 376 /**
 377  * free_points:
 378  * @points: Array of points allocated by parse_points()
 379  *
 380  * Frees all data allocated by parse_points
 381  */
 382 void free_points(GritsPoints *points)
 383 {
 384         gdouble (**_points)[3] = points;
 385         for (int i = 0; _points[i]; i++)
 386                 g_free(_points[i]);
 387         g_free(_points);
 388 }