Merge pull request #1 from jonathonpsmith/add_isg

Introduction of integrated projection handling and ISG grid/geo conversion capability

Author: jonathonpsmith

geodepy/constants.py

@@ -66,6 +66,7 @@ def __init__(self, falseeast, falsenorth, cmscale, zonewidth, initialcm):
 
 
 utm = Projection(500000, 10000000, 0.9996, 6, -177)
+isg = Projection(300000, 5000000, 0.99994, 2, -177)
 
 
 # Helmert 14 parameter transformation

geodepy/convert.py

@@ -8,7 +8,8 @@
 from math import (sin, cos, atan2, radians, degrees,
                   sqrt, cosh, sinh, tan, atan, log)
 import datetime
-from geodepy.constants import utm, grs80
+import warnings
+from geodepy.constants import utm, isg, grs80, ans
 from geodepy.angles import (DECAngle, HPAngle, GONAngle, DMSAngle, DDMAngle,
                             dec2hp, dec2hpa, dec2gon, dec2gona,
                             dec2dms, dec2ddm,
@@ -19,11 +20,6 @@
                             dd2sec, angular_typecheck)
 
 
-# Universal Transverse Mercator Projection Parameters
-# TODO: Add this part into functions: geo2grid, grid2geo, psfandgridconv
-proj = utm
-
-
 def polar2rect(r, theta):
     """
     Converts point in polar coordinates to corresponding rectangular coordinates
@@ -240,7 +236,7 @@ def beta_coeff(ellipsoid):
     return b2, b4, b6, b8, b10, b12, b14, b16
 
 
-def psfandgridconv(xi1, eta1, lat, lon, cm, conf_lat, ellipsoid=grs80):
+def psfandgridconv(xi1, eta1, lat, lon, cm, conf_lat, ellipsoid=grs80, prj=utm):
     """
     Calculates Point Scale Factor and Grid Convergence. Used in convert.geo2grid
     and convert.grid2geo
@@ -270,7 +266,7 @@ def psfandgridconv(xi1, eta1, lat, lon, cm, conf_lat, ellipsoid=grs80):
         p += 2*r * a[r-1] * cos(2*r * xi1) * cosh(2*r * eta1)
         q += 2*r * a[r-1] * sin(2*r * xi1) * sinh(2*r * eta1)
     q = -q
-    psf = (float(proj.cmscale)
+    psf = (float(prj.cmscale)
            * (A / ellipsoid.semimaj)
            * sqrt(q**2 + p**2)
            * ((sqrt(1 + (tan(lat)**2))
@@ -289,7 +285,7 @@ def psfandgridconv(xi1, eta1, lat, lon, cm, conf_lat, ellipsoid=grs80):
     return psf, grid_conv
 
 
-def geo2grid(lat, lon, zone=0, ellipsoid=grs80):
+def geo2grid(lat, lon, zone=0, ellipsoid=grs80, prj=utm):
     """
     Takes a geographic co-ordinate (latitude, longitude) and returns its
     corresponding Hemisphere, Zone and Projection Easting and Northing, Point
@@ -314,23 +310,40 @@ def geo2grid(lat, lon, zone=0, ellipsoid=grs80):
     lon = angular_typecheck(lon)
     # Input Validation - UTM Extents and Values
     zone = int(zone)
-    if zone < 0 or zone > 60:
-        raise ValueError('Invalid Zone - Zones from 1 to 60')
+    if prj == isg:
+        if zone not in (0, 541, 542, 543, 551, 552, 553, 561, 562, 563, 572):
+            raise ValueError('Invalid Zone - Choose from 541, 542, 543, 551, 552, 553, 561, 562, 563, 572')
+    else:
+        if zone < 0 or zone > 60:
+                raise ValueError('Invalid Zone - Zones from 1 to 60')
 
     if lat < -80 or lat > 84:
         raise ValueError('Invalid Latitude - Latitudes from -80 to +84')
 
     if lon < -180 or lon > 180:
         raise ValueError('Invalid Longitude - Longitudes from -180 to +180')
 
+    if prj == isg and ellipsoid != ans:
+        warnings.warn(message='ISG projection should be used with ANS ellipsoid', category=UserWarning)
+
     A = rect_radius(ellipsoid)
     a = alpha_coeff(ellipsoid)
     lat = radians(lat)
     # Calculate Zone
     if zone == 0:
-        zone = int((float(lon) - (
-                    proj.initialcm - (1.5 * proj.zonewidth))) / proj.zonewidth)
-    cm = float(zone * proj.zonewidth) + (proj.initialcm - proj.zonewidth)
+        if prj == isg:
+            amgzone = (float(lon) - (prj.initialcm - (4.5 * prj.zonewidth))) / (prj.zonewidth * 3)
+            subzone = int((amgzone - int(amgzone)) * 3 + 1)
+            amgzone = int(amgzone)
+            zone = int(f'{amgzone}{subzone}')
+        else:
+            zone = int((float(lon) - (prj.initialcm - (1.5 * prj.zonewidth))) / prj.zonewidth)
+    if prj == isg:
+        amgzone = int(str(zone)[:2])
+        subzone = int(str(zone)[2])
+        cm = float((amgzone - 1) * prj.zonewidth * 3 + prj.initialcm + (subzone - 2) * prj.zonewidth)
+    else:
+        cm = float(zone * prj.zonewidth + prj.initialcm - prj.zonewidth)
 
     # Conformal Latitude
     sigx = (ellipsoid.ecc1 * tan(lat)) / sqrt(1 + (tan(lat) ** 2))
@@ -357,14 +370,14 @@ def geo2grid(lat, lon, zone=0, ellipsoid=grs80):
     y = A * xi
 
     # Hemisphere-dependent UTM Projection Co-ordinates
-    east = proj.cmscale * x + proj.falseeast
+    east = prj.cmscale * x + prj.falseeast
     if y < 0:
         hemisphere = 'South'
-        north = proj.cmscale * y + proj.falsenorth
+        north = prj.cmscale * y + prj.falsenorth
     else:
         hemisphere = 'North'
         falsenorth = 0
-        north = proj.cmscale * y + falsenorth
+        north = prj.cmscale * y + falsenorth
 
     # Point Scale Factor and Grid Convergence
     psf, grid_conv = psfandgridconv(xi1, eta1, degrees(lat), lon, cm, conf_lat)
@@ -375,7 +388,7 @@ def geo2grid(lat, lon, zone=0, ellipsoid=grs80):
             round(psf, 8), grid_conv)
 
 
-def grid2geo(zone, east, north, hemisphere='south', ellipsoid=grs80):
+def grid2geo(zone, east, north, hemisphere='south', ellipsoid=grs80, prj=utm):
     """
     Takes a Transverse Mercator grid co-ordinate (Zone, Easting, Northing,
     Hemisphere) and returns its corresponding Geographic Latitude and Longitude,
@@ -393,8 +406,12 @@ def grid2geo(zone, east, north, hemisphere='south', ellipsoid=grs80):
     """
     # Input Validation - UTM Extents and Values
     zone = int(zone)
-    if zone < 0 or zone > 60:
-        raise ValueError('Invalid Zone - Zones from 1 to 60')
+    if prj == isg:
+        if zone not in (541, 542, 543, 551, 552, 553, 561, 562, 563, 572):
+            raise ValueError('Invalid Zone - Choose from 541, 542, 543, 551, 552, 553, 561, 562, 563, 572')
+    else:
+        if zone < 0 or zone > 60:
+            raise ValueError('Invalid Zone - Zones from 1 to 60')
 
     if east < -2830000 or east > 3830000:
         raise ValueError('Invalid Easting - Must be within'
@@ -407,15 +424,18 @@ def grid2geo(zone, east, north, hemisphere='south', ellipsoid=grs80):
     if h != 'north' and h != 'south':
         raise ValueError('Invalid Hemisphere - String, either North or South')
 
+    if prj == isg and ellipsoid != ans:
+        warnings.warn(message='ISG projection should be used with ANS ellipsoid', category=UserWarning)
+
     A = rect_radius(ellipsoid)
     b = beta_coeff(ellipsoid)
     # Transverse Mercator Co-ordinates
-    x = (east - float(proj.falseeast)) / float(proj.cmscale)
+    x = (east - float(prj.falseeast)) / float(prj.cmscale)
     if hemisphere.lower() == 'north':
-        y = -(north / float(proj.cmscale))
+        y = -(north / float(prj.cmscale))
         hemisign = -1
     else:
-        y = (north - float(proj.falsenorth)) / float(proj.cmscale)
+        y = (north - float(prj.falsenorth)) / float(prj.cmscale)
         hemisign = 1
 
     # Transverse Mercator Ratios
@@ -463,7 +483,12 @@ def f1tn(tn, ecc1, ecc1sq):
     lat = degrees(atan(t))
 
     # Compute Longitude
-    cm = float((zone * proj.zonewidth) + proj.initialcm - proj.zonewidth)
+    if prj == isg:
+        amgzone = int(str(zone)[:2])
+        subzone = int(str(zone)[2])
+        cm = float((amgzone - 1) * prj.zonewidth * 3 + prj.initialcm + (subzone - 2) * prj.zonewidth)
+    else:
+        cm = float((zone * prj.zonewidth) + prj.initialcm - prj.zonewidth)
     long_diff = degrees(atan(sinh(eta1) / cos(xi1)))
     long = cm + long_diff
 

geodepy/tests/test_convert.py

@@ -3,6 +3,7 @@
 import datetime
 import numpy as np
 from geodepy.fileio import read_dnacoord
+from geodepy.constants import grs80, ans, isg
 from geodepy.convert import (dec2hp, hp2dec, DMSAngle, DDMAngle, dec2dms,
                              dec2ddm, hp2dms, hp2ddm, dd2sec,
                              yyyydoy_to_date, date_to_yyyydoy, grid2geo,
@@ -80,6 +81,28 @@ def test_geo_grid_transform_interoperability(self):
                                        np.array(test_grid_coords[['east', 'north']].tolist()),
                                        decimal=3)
 
+    def test_geo_grid_transform_interoperability_isg(self):
+        abs_path = os.path.abspath(os.path.dirname(__file__))
+        test_geo_coords = np.genfromtxt(os.path.join(abs_path, 'resources/Test_Conversion_ISG_Geo.csv'),
+                                        delimiter=',',
+                                        dtype='S4,f8,f8',
+                                        names=['site', 'lat', 'lon'])
+
+        test_grid_coords = np.genfromtxt(os.path.join(abs_path, 'resources/Test_Conversion_ISG_Grid.csv'),
+                                         delimiter=',',
+                                         dtype='S4,i4,f8,f8',
+                                         names=['site', 'zone', 'east', 'north'])
+
+        geoed_grid = np.array(list(grid2geo(*x, ellipsoid=ans, prj=isg) for x in test_grid_coords[['zone', 'east', 'north']]))
+        np.testing.assert_almost_equal(geoed_grid[:, :2], np.array(test_geo_coords[['lat', 'lon']].tolist()),
+                                       decimal=8)
+
+        gridded_geo = np.stack(geo2grid(*x, ellipsoid=ans, prj=isg) for x in np.array(test_geo_coords[['lat', 'lon']].tolist()))
+        np.testing.assert_almost_equal(gridded_geo[:, 2:4].astype(float),
+                                       np.array(test_grid_coords[['east', 'north']].tolist()),
+                                       decimal=3)
+
+
     def test_llh2xyz(self):
 
         # Test of single point
@@ -197,6 +220,12 @@ def test_geo2grid(self):
             geo2grid(0, -181, 0)
         with self.assertRaises(ValueError):
             geo2grid(0, 181, 0)
+        # test ValueError raised if not a valid ISG zone
+        with self.assertRaises(ValueError):
+            grid2geo(zone=530, east=300000, north=1500000, ellipsoid=ans, prj=isg)
+        # test UserWarning raised if prj/ellipsoid combination isn't recommended
+        with self.assertWarns(UserWarning):
+            grid2geo(zone=551, east=300000, north=1500000, ellipsoid=grs80, prj=isg)
 
     def test_grid2geo(self):
         abs_path = os.path.abspath(os.path.dirname(__file__))
@@ -233,7 +262,12 @@ def test_grid2geo(self):
             grid2geo(0, 0, 10000001)
         with self.assertRaises(ValueError):
             grid2geo(0, 0, 500000, 'fail')
-
+        # test ValueError raised if not a valid ISG zone
+        with self.assertRaises(ValueError):
+            grid2geo(zone=530, east=300000, north=1500000, ellipsoid=ans, prj=isg)
+        # test UserWarning raised if prj/ellipsoid combination isn't recommended
+        with self.assertWarns(UserWarning):
+            grid2geo(zone=551, east=300000, north=1500000, ellipsoid=grs80, prj=isg)
 
 if __name__ == '__main__':
     unittest.main()