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
75 #include "grits-util.h"
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
87 * Set the latitude, longitude, and elevation for a point.
89 void grits_point_set_lle(GritsPoint *point, gdouble lat, gdouble lon, gdouble elev)
101 * grits_bounds_set_bounds:
107 * Set the north, south, east, and west edges of the bounding box
109 void grits_bounds_set_bounds(GritsBounds *bounds,
110 gdouble n, gdouble s, gdouble e, gdouble w)
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
131 * Convert a point from latitude, longitude, and elevation to x, y and z
134 void lle2xyz(gdouble lat, gdouble lon, gdouble elev,
135 gdouble *x, gdouble *y, gdouble *z)
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);
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
154 * Convert a point from x, y and z coordinates to latitude, longitude, and
157 void xyz2lle(gdouble x, gdouble y, gdouble z,
158 gdouble *lat, gdouble *lon, gdouble *elev)
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);
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
174 * Get the latitude and longitude for a x, y, z value.
176 void xyz2ll(gdouble x, gdouble y, gdouble z,
177 gdouble *lat, gdouble *lon)
179 gdouble rad = sqrt(x*x + y*y + z*z);
180 *lat = incl2lat(acos(y / rad));
181 *lon = azim2lon(atan2(x,z));
186 * @lon_dist: the distance in degrees of longitude
187 * @lat: the latitude to calculate at
189 * Calculate the distance of longitudinal span at a particular latitude.
191 * Returns: the distance in meters
193 gdouble ll2m(gdouble lon_dist, gdouble lat)
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;
203 * @a: the first point
204 * @b: the second point
206 * Calculate the distance between two three dimensional points.
208 * Returns: the distance between the points
210 gdouble distd(gdouble *a, gdouble *b)
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]));
219 * @a: the first longitude
220 * @b: the second longitude
222 * Calculate the average longitude between two longitudes. This is smart about
223 * which side of the globe the resulting longitude is placed on.
225 * Returns: the average
227 gdouble lon_avg(gdouble a, gdouble b)
229 gdouble diff = ABS(a-b);
230 gdouble avg = (a+b)/2;
243 * @b: the center point
244 * @c: the right point
245 * @out: the cross product
247 * Calculate the cross product of three points
249 void crossd3(gdouble *a, gdouble *b, gdouble *c, gdouble *out)
266 * @a: the first vector
267 * @b: the second vector
268 * @c: the cross product
270 * Calculate the cross product of two vectors
272 void crossd(gdouble *a, gdouble *b, gdouble *out)
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];
283 * Calculate the length (magnitude) of a vector.
285 * Returns: the length
287 gdouble lengthd(gdouble *a)
289 return sqrt(a[0] * a[0] +
296 * @a: the first longitude
298 * Normalize a vector to a mangitude of 1
300 void normd(gdouble *a)
302 gdouble total = lengthd(a);
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
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]
323 * parse_points("30,-80 30,-120 50,-120 50,-80", "\t", " ", ",");
325 * Returns: zero-terminated array of groups of points
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)
331 /* Split and count groups */
332 gchar **sgroups = g_strsplit(string, group_sep, -1);
333 int ngroups = g_strv_length(sgroups);
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);
342 /* Create binary coords */
343 gdouble (*coords)[3] = (gpointer)g_new0(gdouble, 3*(ncoords+1));
344 for (int ci = 0; ci < ncoords; ci++) {
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;
359 /* Insert coords into line array */
369 center->lat = (_bounds.n + _bounds.s)/2;
370 center->lon = lon_avg(_bounds.e, _bounds.w);
378 * @points: Array of points allocated by parse_points()
380 * Frees all data allocated by parse_points
382 void free_points(GritsPoints *points)
384 gdouble (**_points)[3] = points;
385 for (int i = 0; _points[i]; i++)