Skip to content

geographic

This module deals with transformation between spherical and Cartesian coordinate systems.

angular_distance(lon1, lat1, lon2, lat2, radians=False)

Compute the great circle distance between two points on the unit sphere. The output is always equivalent to the angle in radians subtended between the vectors connecting the centre of the sphere to the two points.

Note that all input angles are in degrees, unless radians is True.

Parameters:

Name Type Description Default
lon1 float | numpy array

Longitude of first point

 required
lat1 float | numpy array

Latitude of first point

 required
lon2 float | numpy array

Longitude of second point

 required
lat2 float | numpy array

Latitude of second point

 required
radians bool

If True, input angles are in radians; otherwise they are in degrees.

 False

Returns:

Type Description
float | numpy array

angular distance between two points
Source code in terratools/geographic.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def angular_distance(lon1, lat1, lon2, lat2, radians=False):
    """
    Compute the great circle distance between two points on the unit
    sphere.  The output is always equivalent to the angle in radians
    subtended between the vectors connecting the centre of the sphere
    to the two points.

    Note that all input angles are in degrees, unless ``radians`` is
    ``True``.

    :param lon1: Longitude of first point
    :type lon1: float or numpy array
    :param lat1: Latitude of first point
    :type lat1: float or numpy array
    :param lon2: Longitude of second point
    :type lon2: float or numpy array
    :param lat2: Latitude of second point
    :type lat2: float or numpy array
    :param radians: If True, input angles are in radians; otherwise they are
        in degrees.
    :type radians: bool
    :returns: angular distance between two points
    :rtype: float or numpy array
    """
    if not radians:
        lon1 = np.radians(lon1)
        lat1 = np.radians(lat1)
        lon2 = np.radians(lon2)
        lat2 = np.radians(lat2)

    sin_lat1 = np.sin(lat1)
    sin_lat2 = np.sin(lat2)
    cos_lat1 = np.cos(lat1)
    cos_lat2 = np.cos(lat2)
    cos_lon2_lon1 = np.cos(lon2 - lon1)

    distance = np.arctan2(
        np.sqrt(
            (cos_lat2 * np.sin(lon2 - lon1)) ** 2
            + (cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_lon2_lon1) ** 2
        ),
        sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_lon2_lon1,
    )

    return distance

angular_step(lon, lat, azimuth, distance, radians=False)

Compute the final point on the surface of the sphere reached by travelling along a great circle from some starting point along an initial azimuth. The distance covered is in terms of the angle subtended between the start and end points at the centre of the sphere.

Note that all input angles are in degrees, unless radians is True.

Parameters:

Name Type Description Default
lon float | numpy array

Longitude of starting point

 required
lat float | numpy array

Latitude of starting point

 required
azimuth float | numpy array

Azimuth along which to travel from starting point

 required
distance float | numpy array

Great circle distance between starting and final points in terms of angle subtended between them at the centre of the sphere

 required
radians book

If True, input angles are in radians and the output is in radians also; otherwise all input and output are in degrees.

 False

Returns:

Type Description
float, float | two numpy arrays

final longitude and latitude
Source code in terratools/geographic.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def angular_step(lon, lat, azimuth, distance, radians=False):
    """
    Compute the final point on the surface of the sphere reached by travelling
    along a great circle from some starting point along an initial azimuth.
    The distance covered is in terms of the angle subtended between the start
    and end points at the centre of the sphere.

    Note that all input angles are in degrees, unless ``radians`` is
    ``True``.

    :param lon: Longitude of starting point
    :type lon: float or numpy array
    :param lat: Latitude of starting point
    :type lat: float or numpy array
    :param azimuth: Azimuth along which to travel from starting point
    :type azimuth: float or numpy array
    :param distance: Great circle distance between starting and final points
        in terms of angle subtended between them at the centre of the sphere
    :type distance: float or numpy array
    :param radians: If ``True``, input angles are in radians and the output
        is in radians also; otherwise all input and output are in degrees.
    :type radians: book
    :returns: final longitude and latitude
    :rtype: float, float or two numpy arrays
    """
    if not radians:
        lon = np.radians(lon)
        lat = np.radians(lat)
        azimuth = np.radians(azimuth)
        distance = np.radians(distance)

    lat2 = np.arcsin(
        np.sin(lat) * np.cos(distance)
        + np.cos(lat) * np.sin(distance) * np.cos(azimuth)
    )
    lon2 = lon + np.arctan2(
        np.sin(azimuth) * np.sin(distance) * np.cos(lat),
        np.cos(distance) - np.sin(lat) * np.sin(lat2),
    )

    if not radians:
        lon2, lat2 = np.degrees(lon2), np.degrees(lat2)

    return lon2, lat2

azimuth(lon1, lat1, lon2, lat2, radians=False)

Compute the azimuth from point 1 (lon1, lat1) to point 2 (lon2, lat2) on the surface of a sphere.

Note that all input angles are in degrees, unless radians is False.

Parameters:

Name Type Description Default
lon1 float | numpy array

Longitude of first point

 required
lat1 float | numpy array

Latitude of first point

 required
lon2 float | numpy array

Longitude of second point

 required
lat2 float | numpy array

Latitude of second point

 required
radians bool

If True, input angles are in radians; otherwise they are in degrees.

 False

Returns:

Type Description
float | numpy array

angular distance between two points, in radians if radians is True, and in degrees otherwise.
Source code in terratools/geographic.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
def azimuth(lon1, lat1, lon2, lat2, radians=False):
    """
    Compute the azimuth from point 1 (lon1, lat1) to point 2
    (lon2, lat2) on the surface of a sphere.

    Note that all input angles are in degrees, unless ``radians`` is
    ``False``.

    :param lon1: Longitude of first point
    :type lon1: float or numpy array
    :param lat1: Latitude of first point
    :type lat1: float or numpy array
    :param lon2: Longitude of second point
    :type lon2: float or numpy array
    :param lat2: Latitude of second point
    :type lat2: float or numpy array
    :param radians: If True, input angles are in radians; otherwise they are
        in degrees.
    :type radians: bool
    :returns: angular distance between two points, in radians if
        ``radians`` is ``True``, and in degrees otherwise.
    :rtype: float or numpy array
    """
    if not radians:
        lon1 = np.radians(lon1)
        lat1 = np.radians(lat1)
        lon2 = np.radians(lon2)
        lat2 = np.radians(lat2)

    azimuth = np.arctan2(
        np.sin(lon2 - lon1) * np.cos(lat2),
        np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(lon2 - lon1),
    )

    if not radians:
        return np.degrees(azimuth)
    else:
        return azimuth

cart2geog(x, y, z, radians=False)

Convert coordinates in a Cartesian system to a geographic one. If all coordinates are zero, then return all zeros for longitude, latitude and radius.

If all inputs are scalars, output is also scalar; while if all inputs are arrays, then outputs are arrays. Behaviour is undefined if using a mix of scalars and arrays.

The Cartesian system is defined as: - x passes through longitude = 0 and latitude = 0 - y passes through longitude = 90° and latitude = 0 - z passes through the north pole.

Parameters:

Name Type Description Default
x float | numpy array

x coordinate(s)

 required
y float | numpy array

y coordinate(s)

 required
z float | numpy array

z coordinate(s)

 required
radians bool

If True, return output in radians; otherwise return output in degrees.

 False

Returns:

Type Description
3 floats | 3 numpy arrays

longitude, latitude and radius. Longitude and latitude are in degrees, unless radians is True. Radius has the same units as x, y, and z.
Source code in terratools/geographic.py 
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def cart2geog(x, y, z, radians=False):
    """
    Convert coordinates in a Cartesian system to a geographic one.
    If all coordinates are zero, then return all zeros for
    longitude, latitude and radius.

    If all inputs are scalars, output is also scalar; while if all
    inputs are arrays, then outputs are arrays.  Behaviour is
    undefined if using a mix of scalars and arrays.

    The Cartesian system is defined as:
    - x passes through longitude = 0 and latitude = 0
    - y passes through longitude  = 90° and latitude = 0
    - z passes through the north pole.

    :param x: x coordinate(s)
    :type x: float or numpy array
    :param y: y coordinate(s)
    :type y: float or numpy array
    :param z: z coordinate(s)
    :type z: float or numpy array
    :param radians: If ``True``, return output in radians; otherwise
        return output in degrees.
    :type radians: bool
    :returns: longitude, latitude and radius.  Longitude and latitude
        are in degrees, unless ``radians`` is ``True``.  Radius has the
        same units as x, y, and z.
    :rtype: 3 floats or 3 numpy arrays
    """
    r = np.sqrt(x**2 + y**2 + z**2)
    scalar_input = np.isscalar(r)

    if scalar_input and r == 0:
        return 0.0, 0.0, 0.0

    lon = np.arctan2(y, x)
    lat = np.arcsin(z / r)

    # Fix any coordinates where r was zero
    if not scalar_input:
        undef_inds = np.where(r == 0)
        lon[undef_inds] = 0
        lat[undef_inds] = 0
        r[undef_inds] = 0

    if not radians:
        lon = np.degrees(lon)
        lat = np.degrees(lat)

    return lon, lat, r

geog2cart(lon, lat, r, radians=False)

Convert coordinates in a geographic system to a Cartesian system. If all inputs are scalars, output is also scalar; while if all inputs are arrays, then outputs are arrays. Behaviour is undefined if using a mix of scalars and arrays.

The Cartesian system is defined as: - x passes through longitude = 0 and latitude = 0 - y passes through longitude = 90° and latitude = 0 - z passes through the north pole.

Parameters:

Name Type Description Default
lon float | numpy array

Longitude of point(s) in degrees (default), or in radians if radians is True

 required
lat float | numpy array

Latitude of point(s) in degrees (default), or in radians if radians is True

 required
r float | numpy array

Radius of point(s). The units of x, y and z will be the same as those of r.

 required
radians bool

If True, input are in radians; otherwise they are assumed to be in degrees.

 False

Returns:

Type Description
3 floats | 3 numpy arrays

for scalar input, the three coordinates x, y and z; for vector inputs, vectors of x, y and z
Source code in terratools/geographic.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
def geog2cart(lon, lat, r, radians=False):
    """
    Convert coordinates in a geographic system to a Cartesian system.
    If all inputs are scalars, output is also scalar; while if all
    inputs are arrays, then outputs are arrays.  Behaviour is
    undefined if using a mix of scalars and arrays.

    The Cartesian system is defined as:
    - x passes through longitude = 0 and latitude = 0
    - y passes through longitude  = 90° and latitude = 0
    - z passes through the north pole.

    :param lon: Longitude of point(s) in degrees (default), or
        in radians if ``radians`` is ``True``
    :type lon: float or numpy array
    :param lat: Latitude of point(s) in degrees (default), or
        in radians if ``radians`` is ``True``
    :type lat: float or numpy array
    :param r: Radius of point(s).  The units of x, y and z will be
        the same as those of ``r``.
    :type r: float or numpy array
    :param radians: If ``True``, input are in radians; otherwise they
        are assumed to be in degrees.
    :type radians: bool
    :returns: for scalar input, the three coordinates x, y and z; for
        vector inputs, vectors of x, y and z
    :rtype: 3 floats or 3 numpy arrays
    """
    if np.any(r < 0):
        raise ValueError("radius cannot be negative")

    if not radians:
        lon = np.radians(lon)
        lat = np.radians(lat)

    x = r * np.cos(lon) * np.cos(lat)
    y = r * np.sin(lon) * np.cos(lat)
    z = r * np.sin(lat)

    return x, y, z

spherical_triangle_area(lon1, lat1, lon2, lat2, lon3, lat3, r=1, radians=False, tol=None)

Return the spherical area covered by a spherical triangle defined by three geographic coordinates (lon1, lat1), (lon2, lat2) and (lon3, lat3), on a sphere with radius r.

Note: Where two sides of the triangle are similar in length (and the third is about half the other two), l'Huilier's formula becomes unstable. To guard against that case, we use a different expression; see here.

Note that all input angles are in degrees, unless radians is False.

Parameters:

Name Type Description Default
lon1 float | numpy array

Longitude of first point

 required
lat1 float | numpy array

Latitude of first point

 required
lon2 float | numpy array

Longitude of second point

 required
lat2 float | numpy array

Latitude of second point

 required
lon3 float | numpy array

Longitude of third point

 required
lat3 float | numpy array

Latitude of third point

 required
radians bool

Whether input is in radians

 False
tol float

Angle tolerance (always in radians) for collinearity test of three points

 None

Returns:

Type Description
float | numpy array

spherical area of triangle
Source code in terratools/geographic.py 
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
def spherical_triangle_area(
    lon1, lat1, lon2, lat2, lon3, lat3, r=1, radians=False, tol=None
):
    """
    Return the spherical area covered by a spherical triangle defined
    by three geographic coordinates (lon1, lat1), (lon2, lat2) and
    (lon3, lat3), on a sphere with radius r.

    Note: Where two sides of the triangle are similar in length (and the third
    is about half the other two), l'Huilier's formula becomes unstable.  To guard
    against that case, we use a different expression;
    see [here](http://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess).

    Note that all input angles are in degrees, unless ``radians`` is
    ``False``.

    :param lon1: Longitude of first point
    :type lon1: float or numpy array
    :param lat1: Latitude of first point
    :type lat1: float or numpy array
    :param lon2: Longitude of second point
    :type lon2: float or numpy array
    :param lat2: Latitude of second point
    :type lat2: float or numpy array
    :param lon3: Longitude of third point
    :type lon3: float or numpy array
    :param lat3: Latitude of third point
    :type lat3: float or numpy array
    :param radians: Whether input is in radians
    :type radians: bool
    :param tol: Angle tolerance (always in radians) for collinearity test
        of three points
    :type tol: float
    :returns: spherical area of triangle
    :rtype: float or numpy array
    """

    if not radians:
        lon1 = np.radians(lon1)
        lat1 = np.radians(lat1)
        lon2 = np.radians(lon2)
        lat2 = np.radians(lat2)
        lon3 = np.radians(lon3)
        lat3 = np.radians(lat3)

    arc1 = angular_distance(lon2, lat2, lon3, lat3, radians=True)
    arc2 = angular_distance(lon3, lat3, lon1, lat1, radians=True)

    # Angle subtended at point 3
    az1 = azimuth(lon3, lat3, lon1, lat1, radians=True)
    az2 = azimuth(lon3, lat3, lon2, lat2, radians=True)

    # Change in azimuth
    c = np.abs(az2 - az1)
    c = np.minimum(c, 2.0 * np.pi - c)

    tan_arc1 = np.tan(arc1)
    tan_arc2 = np.tan(arc2)

    area = (
        2
        * np.arctan(
            tan_arc1 * tan_arc2 * np.sin(c) / (1 + tan_arc1 * tan_arc2 * np.cos(c))
        )
        * r**2
    )

    return area

triangle_interpolation(lon, lat, lon1, lat1, val1, lon2, lat2, val2, lon3, lat3, val3, radians=False)

Find the interpolated value of a point at (lon, lat), calculated by computing the weighted average of the three surrounding values. Hence the point of interest must lie within the triangle defined by the three other points.

The weighting is computed as follows. The three spherical triangles connecting the point of interest P to the outer triangle's vertices A, B and C are found. Then the area of triangles APB, BPC and CPA are computed, and the weighted average of values a at A, b at B and c at C is computed by giving each value the weight of the area of its opposite triangle, divided by the total area.

            x A (point 1, val1)
           /'\
          / ` \
         /   | \
        /   _+  \
       /__--   ` \
    C x__________x B (point 2, val2)
  (point 3, val3)

Note that all input angles are in degrees, unless radians is False.

Parameters:

Name Type Description Default
lon float | numpy array

Longitude of point of interest within the outer triangle.

 required
lat float | numpy array

Latitude of piont of interest.

 required
lon1 float | numpy array

Longitude of first point

 required
lat1 float | numpy array

Latitude of first point of surrounding triangle

 required
val1 float | numpy array

Value at first point of surrounding triangle

 required
lon2 float | numpy array

Longitude of second point of surrounding triangle

 required
lat2 float | numpy array

Latitude of second point of surrounding triangle

 required
val2 float | numpy array

Value at second point of surrounding triangle

 required
lon3 float | numpy array

Longitude of third point of surrounding triangle

 required
lat3 float | numpy array

Latitude of third point of surrounding triangle

 required
val3 float | numpy array

Value at third point of surrounding triangle

 required
radians bool

Whether input is in radians

 False

Returns:

Type Description
float | numpy array

interpolated value
Source code in terratools/geographic.py 
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
def triangle_interpolation(
    lon, lat, lon1, lat1, val1, lon2, lat2, val2, lon3, lat3, val3, radians=False
):
    """
    Find the interpolated value of a point at (lon, lat), calculated
    by computing the weighted average of the three surrounding values.
    Hence the point of interest must lie within the triangle defined
    by the three other points.

    The weighting is computed as follows.  The three spherical triangles
    connecting the point of interest P to the outer triangle's vertices
    A, B and C are found.  Then the area of triangles APB, BPC and
    CPA are computed, and the weighted average of values a at A, b at B
    and c at C is computed by giving each value the weight of the area
    of its opposite triangle, divided by the total area.

    <pre>
                x A (point 1, val1)
               /'\\
              / ` \\
             /   | \\
            /   _+  \\
           /__--   ` \\
        C x__________x B (point 2, val2)
      (point 3, val3)
    </pre>

    Note that all input angles are in degrees, unless ``radians`` is
    ``False``.

    :param lon: Longitude of point of interest within the outer triangle.
    :type lon: float or numpy array
    :param lat: Latitude of piont of interest.
    :type lat: float or numpy array
    :param lon1: Longitude of first point
    :type lon1: float or numpy array
    :param lat1: Latitude of first point of surrounding triangle
    :type lat1: float or numpy array
    :param val1: Value at first point of surrounding triangle
    :type val1: float or numpy array
    :param lon2: Longitude of second point of surrounding triangle
    :type lon2: float or numpy array
    :param lat2: Latitude of second point of surrounding triangle
    :type lat2: float or numpy array
    :param val2: Value at second point of surrounding triangle
    :type val2: float or numpy array
    :param lon3: Longitude of third point of surrounding triangle
    :type lon3: float or numpy array
    :param lat3: Latitude of third point of surrounding triangle
    :type lat3: float or numpy array
    :param val3: Value at third point of surrounding triangle
    :type val3: float or numpy array
    :param radians: Whether input is in radians
    :type radians: bool
    :returns: interpolated value
    :rtype: float or numpy array
    """
    if not radians:
        lon = np.radians(lon)
        lat = np.radians(lat)
        lon1 = np.radians(lon1)
        lat1 = np.radians(lat1)
        lon2 = np.radians(lon2)
        lat2 = np.radians(lat2)
        lon3 = np.radians(lon3)
        lat3 = np.radians(lat3)

    area12 = spherical_triangle_area(lon, lat, lon1, lat1, lon2, lat2, radians=True)
    area23 = spherical_triangle_area(lon, lat, lon2, lat2, lon3, lat3, radians=True)
    area31 = spherical_triangle_area(lon, lat, lon3, lat3, lon1, lat1, radians=True)

    total_area = area12 + area23 + area31

    interpolated_value = (area23 * val1 + area31 * val2 + area12 * val3) / total_area

    return interpolated_value