Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 7 commits into from
Feb 21, 2022
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c88cb14
define ISG as an instance of a Projection object
jonathonpsmith Feb 10, 2022
740537b
add prj as a parameter to:
jonathonpsmith Feb 10, 2022
5de454e
use the warnings module to warn about mismatch between ISG projection…
jonathonpsmith Feb 16, 2022
4aa33ab
add def test_geo_grid_interoperability_isg() to test ISG functionality
jonathonpsmith Feb 16, 2022
7f99d6b
Merge pull request #1 from jonathonpsmith/add_isg
jonathonpsmith Feb 16, 2022
b100619
test data for ISG<>AGD66 interoperability
jonathonpsmith Feb 18, 2022
d7db0ec
add comments with info and links for ISG
jonathonpsmith Feb 21, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions geodepy/constants.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
73 changes: 49 additions & 24 deletions geodepy/convert.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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
Expand All @@ -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))
Expand All @@ -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)
Expand All @@ -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,
Expand All @@ -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'
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
212 changes: 212 additions & 0 deletions geodepy/tests/resources/Test_Conversion_ISG_Geo.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
Loading