2 * Copyright (C) 2009-2011 Andy Spencer <andy753421@gmail.com>
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.
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.
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/>.
20 * @short_description: Geographic utilities
22 * Miscellaneous utility functions, these deal mostly with coordinate
23 * conversion. Below are some examples that should help demonstrate how these
27 * <title>Terminology</title>
30 * rad - Radians, also radius
31 * m - Meters, for earth-based distances
32 * px - Pixels, for screen-based distances
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
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
43 * x - 0° lon is positive
44 * y - 90° lon is positive
45 * z - North pole is positive
47 * llh - lat,lon,height
56 * <title>Conversions</title>
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
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
74 #include "grits-util.h"
80 * grits_point_set_lle:
81 * @point: the point to modify
82 * @lat: the new latitude
83 * @lon: the new longitude
84 * @elev: the new elevation
86 * Set the latitude, longitude, and elevation for a point.
88 void grits_point_set_lle(GritsPoint *point, gdouble lat, gdouble lon, gdouble elev)
100 * grits_bounds_set_bounds:
106 * Set the north, south, east, and west edges of the bounding box
108 void grits_bounds_set_bounds(GritsBounds *bounds,
109 gdouble n, gdouble s, gdouble e, gdouble w)
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
130 * Convert a point from latitude, longitude, and elevation to x, y and z
133 void lle2xyz(gdouble lat, gdouble lon, gdouble elev,
134 gdouble *x, gdouble *y, gdouble *z)
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);
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
153 * Convert a point from x, y and z coordinates to latitude, longitude, and
156 void xyz2lle(gdouble x, gdouble y, gdouble z,
157 gdouble *lat, gdouble *lon, gdouble *elev)
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);
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
173 * Get the latitude and longitude for a x, y, z value.
175 void xyz2ll(gdouble x, gdouble y, gdouble z,
176 gdouble *lat, gdouble *lon)
178 gdouble rad = sqrt(x*x + y*y + z*z);
179 *lat = incl2lat(acos(y / rad));
180 *lon = azim2lon(atan2(x,z));
185 * @lon_dist: the distance in degrees of longitude
186 * @lat: the latitude to calculate at
188 * Calculate the distance of longitudinal span at a particular latitude.
190 * Returns: the distance in meters
192 gdouble ll2m(gdouble lon_dist, gdouble lat)
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;
202 * @a: the first point
203 * @b: the second point
205 * Calculate the distance between two three dimensional points.
207 * Returns: the distance between the points
209 gdouble distd(gdouble *a, gdouble *b)
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]));
218 * @a: the first longitude
219 * @b: the second longitude
221 * Calculate the average longitude between two longitudes. This is smart about
222 * which side of the globe the resulting longitude is placed on.
224 * Returns: the average
226 gdouble lon_avg(gdouble a, gdouble b)
228 gdouble diff = ABS(a-b);
229 gdouble avg = (a+b)/2;
242 * @b: the center point
243 * @c: the right point
244 * @out: the cross product
246 * Calculate the cross product of three points
248 void crossd3(gdouble *a, gdouble *b, gdouble *c, gdouble *out)
265 * @a: the first vector
266 * @b: the second vector
267 * @c: the cross product
269 * Calculate the cross product of two vectors
271 void crossd(gdouble *a, gdouble *b, gdouble *out)
273 out[0] = a[1] * b[2] - a[2] * b[1];
274 out[1] = a[2] * b[0] - a[0] * b[2];
275 out[2] = a[0] * b[1] - a[1] * b[0];
282 * Calculate the length (magnitude) of a vector.
284 * Returns: the length
286 gdouble lengthd(gdouble *a)
288 return sqrt(a[0] * a[0] +
295 * @a: the first longitude
297 * Normalize a vector to a mangitude of 1
299 void normd(gdouble *a)
301 gdouble total = lengthd(a);