Merge pull request #138 from jonathonpsmith/master

Introduction of a projection parameter for convert.py functions and ISG grid/geo conversion capability

Author: harry093

geodepy/constants.py

@@ -67,6 +67,11 @@ def __init__(self, falseeast, falsenorth, cmscale, zonewidth, initialcm):
 
 utm = Projection(500000, 10000000, 0.9996, 6, -177)
 
+# Integrated Survey Grid - used in NSW as the projection for AGD66
+# Spatial Services projections page - https://www.spatial.nsw.gov.au/surveying/geodesy/projections
+# ISG Technical Manual - https://www.spatial.nsw.gov.au/__data/assets/pdf_file/0017/25730/ISG.pdf
+isg = Projection(300000, 5000000, 0.99994, 2, -177)
+
 
 # Helmert 14 parameter transformation
 class Transformation(object):

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/resources/Test_Conversion_ISG_Geo.csv

@@ -0,0 +1,212 @@
+ALBU,-36.079078850050,146.914301020650,
+ANNA,-32.786387375840,152.085197007870,
+ARD2,-30.522463664480,151.673569431060,
+ARDL,-34.352400821090,146.902223561530,
+ARMD,-30.516201910320,151.664711111400,
+ASHF,-29.325939194990,151.092033399130,
+BALN,-28.874270142210,153.563020333320,
+BANK,-33.916653008040,151.035165581320,
+BARR,-34.565756130890,150.857174374750,
+BATH,-33.431260053750,149.566032786600,
+BBDH,-32.343426762120,146.652773125790,
+BDST,-27.988693001490,152.994004422590,
+BEGA,-36.677477042550,149.840791166900,
+BERM,-36.432192532390,150.074756654240,
+BIGG,-36.282304606530,148.024989032530,
+BING,-32.413362673700,151.651193214890,
+BJCT,-30.103225969580,148.963030107990,
+BKNL,-31.997804118910,141.468816141530,
+BLCK,-31.655869939440,150.243563220700,
+BLKT,-33.773325128960,150.907625503660,
+BLRN,-34.647376511190,143.566122032320,
+BNDC,-37.149042056380,148.883670791220,
+BOMB,-36.913603745640,149.235862495550,
+BOOR,-34.439875740910,148.700493741500,
+BORA,-31.513812898990,150.639630095520,
+BRBA,-30.381927466300,150.606159279950,
+BRDW,-35.448057820140,149.783983323290,
+BRWN,-29.969962049740,146.859629830410,
+BUND,-30.186693592420,151.088547277780,
+BURK,-30.096818183800,145.933090647440,
+CANR,-37.566417946260,149.156254776660,
+CARG,-33.289969788060,146.364698024680,
+CBAR,-31.515200953450,145.834937112410,
+CBLE,-30.955106664800,148.377054074800,
+CBRA,-35.913536753510,145.642887120900,
+CBRM,-35.913730019200,145.644123289160,
+CESS,-32.836839479180,151.355424238120,
+CHCC,-30.296949944770,153.116452056710,
+CHIP,-33.888468395020,151.200058446270,
+CHSM,-35.419814372640,149.119590268640,
+CKWL,-34.457633734950,149.471319973730,
+CLAH,-31.832023492750,149.713778172590,
+CLBI,-29.544125834720,148.584346163510,
+CMDN,-34.047855394520,150.690349966570,
+CNBN,-31.334889140090,149.268408606510,
+CNDA,-30.466878337590,147.686954814500,
+CNDO,-33.086709612700,147.149791898160,
+CNWD,-35.207947226360,149.031553277260,
+COBH,-30.808623541220,146.779606564930,
+COFF,-30.301715186760,153.137240133770,
+COLE,-34.808511810850,145.878982178890,
+COMA,-36.236931036870,149.126002278100,
+COPS,-29.582469259810,152.772422029060,
+CRDX,-34.326381899560,150.766395626120,
+CSNO,-28.867168234440,153.046492333440,
+CTMD,-34.640854520510,148.024428846600,
+CWN2,-33.595280843160,151.170401310160,
+CWRA,-33.832803699280,148.701112901410,
+DBBO,-32.250934874890,148.600941240150,
+DKSN,-35.252310013350,149.134522098020,
+DLQN,-35.533143342320,144.963344877280,
+DORR,-30.351968626410,152.711936850130,
+DUNE,-32.013429920710,149.387007128460,
+DWHY,-33.752851019060,151.285837551170,
+ECHU,-36.141836503430,144.752294114480,
+ECOR,-34.377194893520,150.909653507890,
+EDEN,-37.073313317370,149.908099186570,
+EMMD,-31.667875248990,144.284531118480,
+FLGP,-31.106188458550,141.726864221780,
+FORB,-33.386823485560,148.005764467050,
+FORS,-32.202710872290,152.521279912290,
+FTDN,-33.856641888470,151.224090214610,
+GABO,-37.569690473100,149.913973093290,
+GERR,-34.747858820360,150.826503806750,
+GFEL,-33.894580785140,148.159517594830,
+GFTH,-34.287847702910,146.035544915810,
+GFTN,-29.694628189700,152.931812918630,
+GILG,-31.712549375050,148.661286333990,
+GLB2,-34.747981834810,149.709950157180,
+GLBN,-34.757258766220,149.716521170260,
+GLIN,-29.745140339400,151.750555070350,
+GNGN,-35.186526813120,149.129325278390,
+GNOA,-37.478456864560,149.588255913920,
+GONG,-34.428798726090,150.897673735630,
+GOOL,-33.986013365380,145.705699947480,
+GUNN,-30.979748133500,150.255136195760,
+GURL,-29.736374729300,149.794397024980,
+GWAB,-30.607217796500,148.970055985540,
+HAY1,-34.507065605300,144.851186647770,
+HERN,-30.328217372130,152.485758066980,
+HILL,-33.486454175090,145.530098471910,
+HLB2,-35.712778654860,147.316222432940,
+HLBK,-35.725922870950,147.316155284490,
+HNSB,-33.701901431810,151.096493086520,
+IHOE,-32.865662241970,143.490795121490,
+INVL,-29.778008344970,151.113275953640,
+IRYM,-34.221226895810,142.189912111170,
+IVH2,-32.884967999620,144.307617640270,
+JERI,-35.356857921110,145.724204898170,
+KIRR,-34.044934188350,151.071899350910,
+KNCB,-33.482823083220,151.390789830280,
+KRNG,-35.736914809100,143.921006923710,
+KTMB,-33.697374117700,150.317020257580,
+KULW,-32.327731173920,145.008613436410,
+KURJ,-33.555848060380,150.664437736660,
+LGOW,-33.482521885170,150.158662459500,
+LIPO,-34.106879828990,141.009902661120,
+LIRI,-29.431299550670,147.981730153680,
+LKHT,-35.228801969070,146.704570180350,
+LORD,-31.521522180960,159.060201939740,
+LOTH,-30.535051650830,145.115609446820,
+MACK,-30.712034649050,152.917284371330,
+MENA,-34.127673092760,150.742585967250,
+MENO,-34.273855325240,141.805322115080,
+MGRV,-33.628099764890,150.829824546440,
+MHOP,-32.921796300530,145.910784251350,
+MNDE,-32.394328297570,142.416650264410,
+MNPK,-33.173075032630,151.549959133970,
+MOGO,-35.748958806400,150.190785610250,
+MOUL,-35.092564428440,144.034741161740,
+MREE,-29.459300454610,149.824413821780,
+MRWA,-32.141340244110,150.354104895070,
+MSV2,-34.550684444020,150.372985351130,
+MSVL,-34.552112656670,150.372214490690,
+MTHR,-32.617018805550,151.097961971100,
+MTLD,-32.738569519450,151.559156812120,
+MUDG,-32.591572907010,149.583509963220,
+MWAL,-35.994674877580,145.987190855910,
+NBRI,-30.331701161640,149.785211166790,
+NBRK,-29.678703399780,145.812838240400,
+NDRA,-34.753444164790,146.536290015360,
+NEWE,-32.925559065750,151.787608326110,
+NGAN,-31.565459310400,147.193360053320,
+NMBN,-28.599434985970,153.230499970500,
+NOWE,-31.515142575740,151.715151497200,
+NRMN,-32.236010961490,148.237516634800,
+NSTA,-29.047129217410,150.442979523360,
+NWCS,-32.931168819540,151.764092098380,
+NWRA,-34.875343194560,150.603631645820,
+NYMA,-32.067701391520,146.314409977290,
+OBRN,-33.705521182240,149.856428172010,
+ORNG,-33.286762296710,149.096837893300,
+OVAL,-32.755411057210,148.644863619150,
+OXLY,-34.194693908940,144.100428423170,
+PARK,-33.000328312300,148.263413990430,
+PBOT,-33.975609862890,151.210895108160,
+PCTN,-34.159189725940,150.603976898820,
+PERI,-36.412841838130,148.408734116770,
+PIAN,-35.055844182060,143.326428977960,
+PMAC,-31.463504548720,152.896574077270,
+POON,-33.384043016550,142.565575120360,
+PRCE,-35.365137141950,149.087785486010,
+PRKS,-33.136292183010,148.175265584150,
+PTKL,-34.477127766480,150.912527248830,
+PUTY,-32.954620895850,150.658359913800,
+QUAM,-30.934558432890,147.868825445340,
+RAND,-35.595319266520,146.577251067320,
+RANK,-33.843527001430,146.260407837020,
+RAYM,-32.764467134990,151.744289410730,
+RBVL,-34.591834001070,142.768576347670,
+RGLN,-33.417233215590,149.653804409750,
+ROBI,-28.078561956330,153.380185071210,
+RUUS,-34.043990798410,141.267661676860,
+RYLS,-32.794022742730,149.975508128490,
+SCON,-32.052935807700,150.868499000200,
+SLTC,-28.273328280410,153.573806707950,
+SNGO,-32.559818166330,151.174604604680,
+SPPT,-32.961843433420,151.622519349880,
+SPWD,-33.700119179160,150.562771606160,
+STR1,-35.317094613710,149.008840791660,
+STR2,-35.317722240160,149.008942882160,
+SYDN,-33.782461034370,151.149221129830,
+SYM1,-35.344075493330,149.159851516890,
+TAMW,-31.094437232860,150.929780603990,
+TARE,-31.913796793190,152.462679874520,
+TBOB,-29.451682888950,142.056174284850,
+TID1,-35.400765396740,148.978784924650,
+TLPA,-30.868502342240,144.412558269720,
+TMBA,-35.779238740370,148.010291213980,
+TMRA,-34.448001871050,147.533182888470,
+TMUT,-35.302643610120,148.218953584450,
+TNTR,-29.056362484200,152.018821126120,
+TOTT,-32.254926544980,147.368539197360,
+TULL,-32.633765292060,147.568499137660,
+TURO,-36.036714128190,150.120994043860,
+ULLA,-35.363540535700,150.464154105420,
+UNSW,-33.919253347720,151.230686636390,
+UNW2,-33.919775931790,151.231539899560,
+VLWD,-33.882216612320,150.976004844000,
+WAKL,-35.456861348640,144.379706933050,
+WALW,-35.967558982930,147.732940588930,
+WARI,-29.542163135120,150.573076304270,
+WARW,-28.215068628490,152.029299479150,
+WCHP,-31.462664673680,152.733164012130,
+WDBG,-28.395031300420,152.606047728950,
+WEEM,-29.014858731570,149.253381412470,
+WFAL,-34.135791767390,150.993806632310,
+WGBA,-33.888008907390,150.596577774240,
+WGGA,-35.108730522160,147.368062957900,
+WLCA,-31.556825338460,143.373845848210,
+WLGT,-30.024952808760,148.115771853170,
+WLMD,-35.596106671330,149.136499948870,
+WLWN,-30.932237460500,152.624354982030,
+WRRN,-31.702371232160,147.835254676800,
+WTCF,-30.854815949200,143.091762732270,
+WTON,-32.575358672460,148.955422027630,
+WWLG,-33.704981792290,147.320457917690,
+WYNG,-33.284131475860,151.422832003760,
+YARO,-31.237928944550,151.921103512690,
+YASS,-34.846467014510,148.912040492890,
+YMBA,-29.449038838130,153.356815424290,
+YUNG,-34.305317631760,148.281535497920,