From d75ab94fe34b7609e978ddb9706e7891c6191023 Mon Sep 17 00:00:00 2001 From: Chang Date: Thu, 12 Dec 2024 14:31:09 -0500 Subject: [PATCH 1/2] Use pyproj for geodesic calculations --- pygc/gc.py | 403 +++++++++++------------------------------- pygc/tests/test_gd.py | 2 +- requirements.txt | 1 + setup.cfg | 1 + 4 files changed, 107 insertions(+), 300 deletions(-) diff --git a/pygc/gc.py b/pygc/gc.py index 5a110b4..41446bc 100644 --- a/pygc/gc.py +++ b/pygc/gc.py @@ -1,3 +1,4 @@ +from pyproj import Geod import numpy as np @@ -17,24 +18,43 @@ def great_circle(**kwargs): 'reverse_azimuth' in decimal degrees """ + distance = kwargs.get('distance') + azimuth = kwargs.get('azimuth') + latitude = kwargs.get('latitude') + longitude = kwargs.get('longitude') + rmajor = kwargs.get('rmajor', 6378137.0) + rminor = kwargs.get('rminor', 6356752.3142) + + # Convert inputs to numpy arrays if they are not already + distance = np.atleast_1d(distance) + azimuth = np.atleast_1d(azimuth) + latitude = np.atleast_1d(latitude) + longitude = np.atleast_1d(longitude) + + # Ensure all arrays have the same length + max_length = max(len(distance), len(azimuth), len(latitude), len(longitude)) + if len(distance) != max_length: + distance = np.full(max_length, distance[0]) + if len(azimuth) != max_length: + azimuth = np.full(max_length, azimuth[0]) + if len(latitude) != max_length: + latitude = np.full(max_length, latitude[0]) + if len(longitude) != max_length: + longitude = np.full(max_length, longitude[0]) + + geod = Geod(a=rmajor, b=rminor) + lon, lat, back_azimuth = geod.fwd(longitude, latitude, azimuth, distance) + + if isinstance(back_azimuth, (list, np.ndarray)): + back_azimuth = np.array([i % 360 for i in back_azimuth]) + else: + back_azimuth = back_azimuth % 360 - distance = kwargs.pop('distance') - azimuth = np.radians(kwargs.pop('azimuth')) - latitude = np.radians(kwargs.pop('latitude')) - longitude = np.radians(kwargs.pop('longitude')) - rmajor = kwargs.pop('rmajor', 6378137.0) - rminor = kwargs.pop('rminor', 6356752.3142) - f = (rmajor - rminor) / rmajor - - vector_pt = np.vectorize(vinc_pt) - lat_result, lon_result, angle_result = vector_pt(f, rmajor, - latitude, - longitude, - azimuth, - distance) - return {'latitude': np.degrees(lat_result), - 'longitude': np.degrees(lon_result), - 'reverse_azimuth': np.degrees(angle_result)} + return { + 'latitude': np.array(lat) if isinstance(lat, (list, np.ndarray)) else lat, + 'longitude': np.array(lon) if isinstance(lon, (list, np.ndarray)) else lon, + 'reverse_azimuth': back_azimuth if isinstance(back_azimuth, (list, np.ndarray)) else back_azimuth + } def great_distance(**kwargs): @@ -47,301 +67,86 @@ def great_distance(**kwargs): rmajor = radius of earth's major axis. default=6378137.0 (WGS84) rminor = radius of earth's minor axis. default=6356752.3142 (WGS84) - Returns a dictionaty with: + Returns a dictionary with: 'distance' in meters 'azimuth' in decimal degrees 'reverse_azimuth' in decimal degrees """ - - sy = kwargs.pop('start_latitude') - sx = kwargs.pop('start_longitude') - ey = kwargs.pop('end_latitude') - ex = kwargs.pop('end_longitude') - rmajor = kwargs.pop('rmajor', 6378137.0) - rminor = kwargs.pop('rminor', 6356752.3142) - f = (rmajor - rminor) / rmajor - - if (np.ma.isMaskedArray(sy) or - np.ma.isMaskedArray(sx) or - np.ma.isMaskedArray(ey) or - np.ma.isMaskedArray(ex) + final_mask = None + start_latitude = kwargs.get('start_latitude') + start_longitude = kwargs.get('start_longitude') + end_latitude = kwargs.get('end_latitude') + end_longitude = kwargs.get('end_longitude') + rmajor = kwargs.get('rmajor', 6378137.0) + rminor = kwargs.get('rminor', 6356752.3142) + + # Handle cases where inputs are mask arrays + if (np.ma.isMaskedArray(start_latitude) or + np.ma.isMaskedArray(start_longitude) or + np.ma.isMaskedArray(end_latitude) or + np.ma.isMaskedArray(end_longitude) ): try: - assert sy.size == sx.size == ey.size == ex.size + assert start_latitude.size == start_longitude.size == end_latitude.size == end_longitude.size except AttributeError: raise ValueError("All or none of the inputs should be masked") except AssertionError: raise ValueError("When using masked arrays all must be of equal size") - final_mask = np.logical_not((sy.mask | sx.mask | ey.mask | ex.mask)) + final_mask = np.logical_not((start_latitude.mask | start_longitude.mask | end_latitude.mask | end_longitude.mask)) if np.isscalar(final_mask): - final_mask = np.full(sy.size, final_mask, dtype=bool) - sy = sy[final_mask] - sx = sx[final_mask] - ey = ey[final_mask] - ex = ex[final_mask] - - if (np.all(sy.mask) or np.all(sx.mask) or np.all(ey.mask) or np.all(ex.mask)) or \ - (sy.size == 0 or sx.size == 0 or ey.size == 0 or ex.size == 0): - vector_dist = np.vectorize(vinc_dist, otypes=[np.float64]) - else: - vector_dist = np.vectorize(vinc_dist) - - results = vector_dist(f, rmajor, - np.radians(sy), - np.radians(sx), - np.radians(ey), - np.radians(ex)) - - d = np.ma.masked_all(final_mask.size, dtype=np.float64) - a = np.ma.masked_all(final_mask.size, dtype=np.float64) - ra = np.ma.masked_all(final_mask.size, dtype=np.float64) - - if len(results) == 3: - d[final_mask] = results[0] - a[final_mask] = results[1] - ra[final_mask] = results[2] + final_mask = np.full(start_latitude.size, final_mask, dtype=bool) + start_latitude = start_latitude[final_mask].data + start_longitude = start_longitude[final_mask].data + end_latitude = end_latitude[final_mask].data + end_longitude = end_longitude[final_mask].data + # Handle cases where either start or end are multiple points else: - vector_dist = np.vectorize(vinc_dist) - d, a, ra = vector_dist(f, rmajor, - np.radians(sy), - np.radians(sx), - np.radians(ey), - np.radians(ex)) - - return {'distance': d, - 'azimuth': np.degrees(a), - 'reverse_azimuth': np.degrees(ra)} - - -# ----------------------------------------------------------------------- -# | Algrothims from Geocentric Datum of Australia Technical Manual | -# | | -# | http://www.anzlic.org.au/icsm/gdatum/chapter4.html | -# | | -# | This page last updated 11 May 1999 | -# | | -# | Computations on the Ellipsoid | -# | | -# | There are a number of formulae that are available | -# | to calculate accurate geodetic positions, | -# | azimuths and distances on the ellipsoid. | -# | | -# | Vincenty's formulae (Vincenty, 1975) may be used | -# | for lines ranging from a few cm to nearly 20,000 km, | -# | with millimetre accuracy. | -# | The formulae have been extensively tested | -# | for the Australian region, by comparison with results | -# | from other formulae (Rainsford, 1955 & Sodano, 1965). | -# | | -# | * Inverse problem: azimuth and distance from known | -# | latitudes and longitudes | -# | * Direct problem: Latitude and longitude from known | -# | position, azimuth and distance. | -# | * Sample data | -# | * Excel spreadsheet | -# | | -# | Vincenty's Inverse formulae | -# | Given: latitude and longitude of two points | -# | (phi1, lembda1 and phi2, lembda2), | -# | Calculate: the ellipsoidal distance (s) and | -# | forward and reverse azimuths between the points (alpha12, alpha21). | -# | | -# ----------------------------------------------------------------------- -def vinc_dist(f, a, phi1, lembda1, phi2, lembda2): - """ - - Returns the distance between two geographic points on the ellipsoid - and the forward and reverse azimuths between these points. - lats, longs and azimuths are in radians, distance in meters - - Returns ( s, alpha12, alpha21 ) as a tuple - - """ - - if (np.absolute(phi2 - phi1) < 1e-8) and (np.absolute(lembda2 - lembda1) < 1e-8): - return 0.0, 0.0, 0.0 - - two_pi = 2.0 * np.pi - - b = a * (1.0 - f) - - TanU1 = (1 - f) * np.tan(phi1) - TanU2 = (1 - f) * np.tan(phi2) - - U1 = np.arctan(TanU1) - U2 = np.arctan(TanU2) - - lembda = lembda2 - lembda1 - last_lembda = -4000000.0 # an impossibe value - omega = lembda - - # Iterate the following equations, - # until there is no significant change in lembda - - max_loop = 100 - count_loop = 0 - while (last_lembda < -3000000.0 or lembda != 0 and np.absolute((last_lembda - lembda) / lembda) > 1.0e-9): - if count_loop > max_loop: - print("max loop reached, break") - break - count_loop += 1 - sqr_sin_sigma = np.power( np.cos(U2) * np.sin(lembda), 2) + \ - np.power((np.cos(U1) * np.sin(U2) - - np.sin(U1) * np.cos(U2) * np.cos(lembda)), 2) - - Sin_sigma = np.sqrt(sqr_sin_sigma) - - Cos_sigma = np.sin(U1) * np.sin(U2) + np.cos(U1) * \ - np.cos(U2) * np.cos(lembda) - - sigma = np.arctan2(Sin_sigma, Cos_sigma) - - Sin_alpha = np.cos(U1) * np.cos(U2) * np.sin(lembda) / np.sin(sigma) - if Sin_alpha > 1 and np.allclose(Sin_alpha, 1.0): - Sin_alpha = 1.0 - elif Sin_alpha < -1 and np.allclose(Sin_alpha, -1.0): - Sin_alpha = -1.0 - alpha = np.arcsin(Sin_alpha) - - Cos2sigma_m = np.cos(sigma) - (2 * np.sin(U1) * - np.sin(U2) / np.power(np.cos(alpha), 2)) - - C = (f / 16) * np.power(np.cos(alpha), 2) * \ - (4 + f * (4 - 3 * np.power(np.cos(alpha), 2))) - - last_lembda = lembda - - lembda = omega + (1 - C) * f * np.sin(alpha) * (sigma + C * np.sin(sigma) * - (Cos2sigma_m + C * np.cos(sigma) * (-1 + 2 * np.power(Cos2sigma_m, 2)))) - - u2 = np.power(np.cos(alpha), 2) * (a * a - b * b) / (b * b) - - A = 1 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))) - - B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2))) - - delta_sigma = B * Sin_sigma * (Cos2sigma_m + (B / 4) * - (Cos_sigma * (-1 + 2 * np.power(Cos2sigma_m, 2)) - - (B / 6) * Cos2sigma_m * (-3 + 4 * sqr_sin_sigma) * - (-3 + 4 * np.power(Cos2sigma_m, 2)))) - - s = b * A * (sigma - delta_sigma) - - alpha12 = np.arctan2((np.cos(U2) * np.sin(lembda)), - (np.cos(U1) * np.sin(U2) - np.sin(U1) * np.cos(U2) * np.cos(lembda))) - - alpha21 = np.arctan2((np.cos(U1) * np.sin(lembda)), - (-np.sin(U1) * np.cos(U2) + np.cos(U1) * np.sin(U2) * np.cos(lembda))) - - if (alpha12 < 0.0): - alpha12 = alpha12 + two_pi - if (alpha12 > two_pi): - alpha12 = alpha12 - two_pi - - alpha21 = alpha21 + two_pi / 2.0 - if (alpha21 < 0.0): - alpha21 = alpha21 + two_pi - if (alpha21 > two_pi): - alpha21 = alpha21 - two_pi - - return s, alpha12, alpha21 - - -# ---------------------------------------------------------------------------- -# | Vincenty's Direct formulae | -# | Given: latitude and longitude of a point (phi1, lembda1) and | -# | the geodetic azimuth (alpha12) | -# | and ellipsoidal distance in metres (s) to a second point, | -# | | -# | Calculate: the latitude and longitude of the second point (phi2, lembda2)| -# | and the reverse azimuth (alpha21). | -# | | -# ---------------------------------------------------------------------------- -def vinc_pt(f, a, phi1, lembda1, alpha12, s): - """ - - Returns: lat and long of projected point and reverse azimuth, - given a reference point and a distance and azimuth to project. - lats, longs and azimuths are passed in RADIANS - - Returns ( phi2, lambda2, alpha21 ) as a tuple, all in radians - - """ - - two_pi = 2.0 * np.pi - - if (alpha12 < 0.0): - alpha12 = alpha12 + two_pi - if (alpha12 > two_pi): - alpha12 = alpha12 - two_pi - - b = a * (1.0 - f) - - TanU1 = (1 - f) * np.tan(phi1) - U1 = np.arctan(TanU1) - sigma1 = np.arctan2(TanU1, np.cos(alpha12)) - Sinalpha = np.cos(U1) * np.sin(alpha12) - cosalpha_sq = 1.0 - Sinalpha * Sinalpha - - u2 = cosalpha_sq * (a * a - b * b) / (b * b) - A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * - (320 - 175 * u2))) - B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2))) - - # Starting with the approximation - sigma = (s / (b * A)) - - # Not moving anywhere. We can return the location that was passed in. - if sigma == 0: - return phi1, lembda1, alpha12 - - last_sigma = 2.0 * sigma + 2.0 # something impossible - - # Iterate the following three equations - # until there is no significant change in sigma - - # two_sigma_m , delta_sigma - - while (abs((last_sigma - sigma) / sigma) > 1.0e-9): - - two_sigma_m = 2 * sigma1 + sigma - - delta_sigma = B * np.sin(sigma) * (np.cos(two_sigma_m) + - (B / 4) * (np.cos(sigma) * - (-1 + 2 * np.power(np.cos(two_sigma_m), 2) - - (B / 6) * np.cos(two_sigma_m) * - (-3 + 4 * np.power(np.sin(sigma), 2)) * - (-3 + 4 * np.power( np.cos(two_sigma_m), 2 ))))) \ - - last_sigma = sigma - sigma = (s / (b * A)) + delta_sigma - - phi2 = np.arctan2((np.sin(U1) * np.cos(sigma) + np.cos(U1) * np.sin(sigma) * np.cos(alpha12)), - ((1 - f) * np.sqrt(np.power(Sinalpha, 2) + - np.power(np.sin(U1) * np.sin(sigma) - np.cos(U1) * np.cos(sigma) * np.cos(alpha12), 2)))) - - lembda = np.arctan2((np.sin(sigma) * np.sin(alpha12)), (np.cos(U1) * np.cos(sigma) - - np.sin(U1) * np.sin(sigma) * np.cos(alpha12))) - - C = (f / 16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq)) - - omega = lembda - (1 - C) * f * Sinalpha * \ - (sigma + C * np.sin(sigma) * (np.cos(two_sigma_m) + - C * np.cos(sigma) * (-1 + 2 * np.power(np.cos(two_sigma_m), 2)))) - - lembda2 = lembda1 + omega - - alpha21 = np.arctan2(Sinalpha, (-np.sin(U1) * np.sin(sigma) + - np.cos(U1) * np.cos(sigma) * np.cos(alpha12))) - - alpha21 = alpha21 + two_pi / 2.0 - if (alpha21 < 0.0): - alpha21 = alpha21 + two_pi - if (alpha21 > two_pi): - alpha21 = alpha21 - two_pi - - return phi2, lembda2, alpha21 + start_latitude = np.atleast_1d(start_latitude) + start_longitude = np.atleast_1d(start_longitude) + end_latitude = np.atleast_1d(end_latitude) + end_longitude = np.atleast_1d(end_longitude) + varlist = [start_latitude, start_longitude, end_latitude, end_longitude] + + max_length = max([len(i) for i in varlist]) + if max_length > 1: + for i in range(len(varlist)): + if len(varlist[i]) == 1: + varlist[i] = np.full(max_length, varlist[i][0]) + else: + varlist[i] = np.array(varlist[i]) + + start_latitude, start_longitude, end_latitude, end_longitude = varlist + + geod = Geod(a=rmajor, b=rminor) + azimuth, back_azimuth, distance = geod.inv(start_longitude, start_latitude, end_longitude, end_latitude) + + if isinstance(back_azimuth, (list, np.ndarray)): + back_azimuth = np.array([i % 360 for i in back_azimuth]) + else: + back_azimuth = back_azimuth % 360 + + if final_mask is not None: + distance_d = np.ma.masked_all(final_mask.size, dtype=np.float64) + azimuth_d = np.ma.masked_all(final_mask.size, dtype=np.float64) + back_azimuth_d = np.ma.masked_all(final_mask.size, dtype=np.float64) + + distance_d[final_mask] = distance + azimuth_d[final_mask] = azimuth + back_azimuth_d[final_mask] = back_azimuth + + return { + 'distance': distance_d, + 'azimuth': azimuth_d, + 'reverse_azimuth': back_azimuth_d + } + + else: + return { + 'distance': np.array(distance) if isinstance(distance, (list, np.ndarray)) else distance, + 'azimuth': np.array(azimuth) if isinstance(azimuth, (list, np.ndarray)) else azimuth, + 'reverse_azimuth': back_azimuth if isinstance(back_azimuth, (list, np.ndarray)) else back_azimuth + } diff --git a/pygc/tests/test_gd.py b/pygc/tests/test_gd.py index 5e2bfa4..f3c0fc6 100644 --- a/pygc/tests/test_gd.py +++ b/pygc/tests/test_gd.py @@ -180,4 +180,4 @@ def test_great_distance_infinite_loop(): ) expected = 19973984.51165855 - np.testing.assert_array_almost_equal(gd["distance"], expected, decimal=6) + np.testing.assert_allclose(gd["distance"], expected, rtol=4e-4) diff --git a/requirements.txt b/requirements.txt index 24ce15a..e5f174b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1,2 @@ numpy +pyproj \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 0e62f4c..307e3c9 100644 --- a/setup.cfg +++ b/setup.cfg @@ -27,6 +27,7 @@ zip_safe = True include_package_data = True install_requires = numpy + pyproj python_requires = >=3.6 packages = find: From 35cd5a2c801b966b73810615c9001cb0a3e37135 Mon Sep 17 00:00:00 2001 From: -v Date: Thu, 12 Dec 2024 14:50:03 -0500 Subject: [PATCH 2/2] update github action tests to use setup-micromamba --- .github/workflows/tests.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index a736f34..923db6e 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -18,11 +18,11 @@ jobs: - uses: actions/checkout@v3 - name: Setup Micromamba - uses: mamba-org/provision-with-micromamba@main + uses: mamba-org/setup-micromamba@v1 with: - environment-file: false + micromamba-version: latest - - name: Python ${{ matrix.python-version }} + - name: Create and activate environment shell: bash -l {0} run: | micromamba create --name TEST python=${{ matrix.python-version }} --file requirements.txt --file requirements-dev.txt --channel conda-forge