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  * GisPoint *
  78  ************/
  79 /**
  80  * gis_point_set_lle:
  81  * @point: the point to modify
  82  * @lat:   the new latitude
  83  * @lon:   the new longitude
  84  * @elev:  the new elevation
  85  *
  86  * Set the latitude, longitude, and elevation for a point.
  87  */
  88 void gis_point_set_lle(GisPoint *point, gdouble lat, gdouble lon, gdouble elev)
  89 {
  90         point->lat  = lat;
  91         point->lon  = lon;
  92         point->elev = elev;
  93 }
  94
  95
  96 /***********
  97  * GisBBox *
  98  ***********/
  99 /**
 100  * gis_bbox_set_bounds:
 101  * @n: the north edge
 102  * @s: the south edge
 103  * @e: the east edge
 104  * @w: the west edge
 105  *
 106  * Set the north, south, east, and west edges of the bounding box
 107  */
 108 void gis_bbox_set_bounds(GisBBox *bbox,
 109                 gdouble n, gdouble s, gdouble e, gdouble w)
 110 {
 111         bbox->n = n;
 112         bbox->s = s;
 113         bbox->e = e;
 114         bbox->w = w;
 115 }
 116
 117
 118 /******************
 119  * Global helpers *
 120  ******************/
 121 /**
 122  * lle2xyz:
 123  * @lat:  the latitude
 124  * @lon:  the longitude
 125  * @elev: the elevation
 126  * @x:    the resulting x coordinate
 127  * @y:    the resulting y coordinate
 128  * @z:    the resulting z coordinate
 129  *
 130  * Convert a point from latitude, longitude, and elevation to x, y and z
 131  * coordinates.
 132  */
 133 void lle2xyz(gdouble lat, gdouble lon, gdouble elev,
 134                 gdouble *x, gdouble *y, gdouble *z)
 135 {
 136         gdouble rad  = elev2rad(elev);
 137         gdouble azim = lon2azim(lon);
 138         gdouble incl = lat2incl(lat);
 139         *z = rad * cos(azim) * sin(incl);
 140         *x = rad * sin(azim) * sin(incl);
 141         *y = rad * cos(incl);
 142 }
 143
 144 /**
 145  * xyz2lle:
 146  * @x:    the x coordinate
 147  * @y:    the y coordinate
 148  * @z:    the z coordinate
 149  * @lat:  the resulting latitude
 150  * @lon:  the resulting longitude
 151  * @elev: the resulting elevation
 152  *
 153  * Convert a point from x, y and z coordinates to latitude, longitude, and
 154  * elevation.
 155  */
 156 void xyz2lle(gdouble x, gdouble y, gdouble z,
 157                 gdouble *lat, gdouble *lon, gdouble *elev)
 158 {
 159         gdouble rad = sqrt(x*x + y*y + z*z);
 160         *lat  = incl2lat(acos(y / rad));
 161         *lon  = azim2lon(atan2(x,z));
 162         *elev = rad2elev(rad);
 163 }
 164
 165 /**
 166  * xyz2ll:
 167  * @x:    the x coordinate
 168  * @y:    the y coordinate
 169  * @z:    the z coordinate
 170  * @lat:  the resulting latitude
 171  * @lon:  the resulting longitude
 172  *
 173  * Get the latitude and longitude for a x, y, z value.
 174  */
 175 void xyz2ll(gdouble x, gdouble y, gdouble z,
 176                 gdouble *lat, gdouble *lon)
 177 {
 178         gdouble rad = sqrt(x*x + y*y + z*z);
 179         *lat = incl2lat(acos(y / rad));
 180         *lon = azim2lon(atan2(x,z));
 181 }
 182
 183 /**
 184  * ll2m:
 185  * @lon_dist: the distance in degrees of longitude
 186  * @lat:      the latitude to calculate at
 187  *
 188  * Calculate the distance of longitudinal span at a particular latitude.
 189  *
 190  * Returns: the distance in meters
 191  */
 192 gdouble ll2m(gdouble lon_dist, gdouble lat)
 193 {
 194         gdouble incl = lat2incl(lat);
 195         gdouble rad  = sin(incl) * EARTH_R;
 196         gdouble circ = 2 * G_PI * rad;
 197         return lon_dist/360 * circ;
 198 }
 199
 200 /**
 201  * distd:
 202  * @a: the first point
 203  * @b: the second point
 204  *
 205  * Calculate the distance between two three dimensional points.
 206  *
 207  * Returns: the distance between the points
 208  */
 209 gdouble distd(gdouble *a, gdouble *b)
 210 {
 211         return sqrt((a[0]-b[0])*(a[0]-b[0]) +
 212                     (a[1]-b[1])*(a[1]-b[1]) +
 213                     (a[2]-b[2])*(a[2]-b[2]));
 214 }
 215
 216 /**
 217  * lon_avg:
 218  * @a: the first longitude
 219  * @b: the second longitude
 220  *
 221  * Calculate the average longitude between two longitudes. This is smart about
 222  * which side of the globe the resulting longitude is placed on.
 223  *
 224  * Returns: the average
 225  */
 226 gdouble lon_avg(gdouble a, gdouble b)
 227 {
 228         gdouble diff = ABS(a-b);
 229         gdouble avg  = (a+b)/2;
 230         if (diff > 180) {
 231                 if (avg >= 0)
 232                         avg -= 180;
 233                 else
 234                         avg += 180;
 235         }
 236         return avg;
 237 }