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

   1 /*
   2  * Copyright (C) 2009-2010 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:gis-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
  74 #include "gis-util.h"
  75
  76 /******************
  77  * Global helpers *
  78  ******************/
  79 /**
  80  * lle2xyz:
  81  * @lat:  the latitude
  82  * @lon:  the longitude
  83  * @elev: the elevation
  84  * @x:    the resulting x coordinate
  85  * @y:    the resulting y coordinate
  86  * @z:    the resulting z coordinate
  87  *
  88  * Convert a point from latitude, longitude, and elevation to x, y and z
  89  * coordinates.
  90  */
  91 void lle2xyz(gdouble lat, gdouble lon, gdouble elev,
  92                 gdouble *x, gdouble *y, gdouble *z)
  93 {
  94         gdouble rad  = elev2rad(elev);
  95         gdouble azim = lon2azim(lon);
  96         gdouble incl = lat2incl(lat);
  97         *z = rad * cos(azim) * sin(incl);
  98         *x = rad * sin(azim) * sin(incl);
  99         *y = rad * cos(incl);
 100 }
 101
 102 /**
 103  * xyz2lle:
 104  * @x:    the x coordinate
 105  * @y:    the y coordinate
 106  * @z:    the z coordinate
 107  * @lat:  the resulting latitude
 108  * @lon:  the resulting longitude
 109  * @elev: the resulting elevation
 110  *
 111  * Convert a point from x, y and z coordinates to latitude, longitude, and
 112  * elevation.
 113  */
 114 void xyz2lle(gdouble x, gdouble y, gdouble z,
 115                 gdouble *lat, gdouble *lon, gdouble *elev)
 116 {
 117         gdouble rad = sqrt(x*x + y*y + z*z);
 118         *lat  = incl2lat(acos(y / rad));
 119         *lon  = azim2lon(atan2(x,z));
 120         *elev = rad2elev(rad);
 121 }
 122
 123 /**
 124  * xyz2ll:
 125  * @x:    the x coordinate
 126  * @y:    the y coordinate
 127  * @z:    the z coordinate
 128  * @lat:  the resulting latitude
 129  * @lon:  the resulting longitude
 130  *
 131  * Get the latitude and longitude for a x, y, z value.
 132  */
 133 void xyz2ll(gdouble x, gdouble y, gdouble z,
 134                 gdouble *lat, gdouble *lon)
 135 {
 136         gdouble rad = sqrt(x*x + y*y + z*z);
 137         *lat = incl2lat(acos(y / rad));
 138         *lon = azim2lon(atan2(x,z));
 139 }
 140
 141 /**
 142  * ll2m:
 143  * @lon_dist: the distance in degrees of longitude
 144  * @lat:      the latitude to calculate at
 145  *
 146  * Calculate the distance of longitudinal span at a particular latitude.
 147  *
 148  * Returns: the distance in meters
 149  */
 150 gdouble ll2m(gdouble lon_dist, gdouble lat)
 151 {
 152         gdouble azim = (-lat+90)/180*M_PI;
 153         gdouble rad  = sin(azim) * EARTH_R;
 154         gdouble circ = 2 * M_PI * rad;
 155         return lon_dist/360 * circ;
 156 }
 157
 158 /**
 159  * distd:
 160  * @a: the first point
 161  * @b: the second point
 162  *
 163  * Calculate the distance between two three dimensional points.
 164  *
 165  * Returns: the distance between the points
 166  */
 167 gdouble distd(gdouble *a, gdouble *b)
 168 {
 169         return sqrt((a[0]-b[0])*(a[0]-b[0]) +
 170                     (a[1]-b[1])*(a[1]-b[1]) +
 171                     (a[2]-b[2])*(a[2]-b[2]));
 172 }
 173
 174 /**
 175  * lon_avg:
 176  * @a: the first longitude
 177  * @b: the second longitude
 178  *
 179  * Calculate the average longitude between two longitudes. This is smart about
 180  * which side of the globe the resulting longitude is placed on.
 181  *
 182  * Returns: the average
 183  */
 184 gdouble lon_avg(gdouble a, gdouble b)
 185 {
 186         gdouble diff = ABS(a-b);
 187         gdouble avg  = (a+b)/2;
 188         if (diff > 180) {
 189                 if (avg >= 0)
 190                         avg -= 180;
 191                 else
 192                         avg += 180;
 193         }
 194         return avg;
 195 }