Adding Height Transformation code

.travis.yml

@@ -22,7 +22,7 @@ install:
   # Useful for debugging any issues with conda
   - conda info -a
 
-  - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy
+  - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy gdal
   - source activate test-environment
   - pip install pytest pytest-cov coveralls
 

geodepy/Height_filenames.py

@@ -0,0 +1,12 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 29 08:56:15 2019
+
+@author: u84157
+"""
+file_DOV_PV='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/DOV_PV.tif'
+file_DOV_PM='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/DOV_PM.tif'
+file_AG2020='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AUSGeoid2020_RELEASEV20170908.tif'
+file_AG2020_STD='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AUSGeoid2020_RELEASEV20170908_err.tif'
+file_AVWS='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AVWS_20191107.tif'
+file_AVWS_STD='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AVWS_STD_20191107.tif'

geodepy/heights_s3.py

@@ -0,0 +1,90 @@
+#___________________________________________________________________________#
+# Some notes:
+# Written by Jack McCubbine of Geoscience Australia, date: 08/11/2019
+# This code contains functions to handle tranformations between GPS and 
+# AWVS/AHD and Vice Versa 
+# Gridded data used for the varisous reference surfaces are geotiff files
+# These allow direct access remotely using "gdal"
+#___________________________________________________________________________#
+# Import dependencies 
+import gdal
+import numpy as np
+from scipy.interpolate import griddata
+import geodepy.Height_filenames as Height_filenames
+#___________________________________________________________________________#
+# Interpolation functions
+def interp_file(Lat,Long,file):
+    # Import the DOVPM file
+    f = gdal.Open(file)
+    # load band (akin to a variable in dataset)
+    band = f.GetRasterBand(1)   
+    # get the pixel width, height, etc.
+    transform = f.GetGeoTransform()
+    # Grid resolution (known)
+    res=transform[1]
+    # convert lat,lon to row,col
+    column = (Long - transform[0]) / transform[1]
+    row = (Lat - transform[3]) / transform[5]
+    # get pixel values surrounding data point
+    Surrounding_data=(band.ReadAsArray(np.floor(column-2), np.floor(row-2), 5, 5))
+    # convert row,col back to north,east
+    Long_c = transform[0] + np.floor(column) * res
+    Lat_c = transform[3] - np.floor(row) * res
+    # set up matrices for interpolation
+    count=-1
+    pos=np.zeros((25,2))
+    Surrounding_data_v=np.zeros((25,1))
+    for k in range(-2,3):
+        for j in range(-2,3):
+            count=count+1
+            pos[count]=(Long_c+j*res,Lat_c-k*res)
+            Surrounding_data_v[count]=Surrounding_data[k+2,j+2]          
+    interp_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic')
+    return interp_val
+
+#___________________________________________________________________________#
+# Functions to handle the conversions from one height to another
+def GPS_to_AVWS(Lat,Long,GPS_H):
+    zeta=interp_file(Lat,Long,Height_filenames.file_AVWS) # AVWS file
+    zeta_std=interp_file(Lat,Long,Height_filenames.file_AVWS_STD) # AVWS STD file
+    NORMAL_H=GPS_H-zeta
+    return [NORMAL_H,zeta_std]
+
+def AVWS_to_GPS(Lat,Long,AVWS_H):
+    zeta=interp_file(Lat,Long,Height_filenames.file_AVWS) # AVWS file
+    zeta_std=interp_file(Lat,Long,Height_filenames.file_AVWS_STD) # AVWS STD file
+    GPS_H=AVWS_H+zeta
+    return [GPS_H,zeta_std]
+
+def AHD_to_AVWS(Lat,Long,AHD_H,file_AG2020):
+    # Convert to GPS
+    GPS_H=AHD_H+interp_file(Lat,Long,Height_filenames.file_AG2020) # AUSGEOID2020 file
+    # Convert to AVWS
+    Normal_H=GPS_H-interp_file(Lat,Long,Height_filenames.file_AVWS) # AVWS file
+    return [Normal_H]
+
+def GPS_to_AHD(Lat,Long,GPS_H):
+    N=interp_file(Lat,Long,Height_filenames.file_AG2020) # AUSGEOID2020 file
+    N_std=interp_file(Lat,Long,Height_filenames.file_AG2020_STD) # AUSGEOID2020 STD file
+    AHD_H=GPS_H-N
+    return [AHD_H,N_std]
+
+def AHD_to_GPS(Lat,Long,AHD_H):
+    N=interp_file(Lat,Long,Height_filenames.file_AG2020) # AUSGEOID2020 file
+    N_std=interp_file(Lat,Long,Height_filenames.file_AG2020_STD) # AUSGEOID2020 STD file
+    GPS_H=AHD_H+N
+    return [GPS_H,N_std]
+
+def AVWS_to_AHD(Lat,Long,Normal_H):   
+    # Convert to GPS    
+    GPS_H=Normal_H+interp_file(Lat,Long,Height_filenames.file_AVWS) # AVWS file
+    # Convert to AHD
+    AHD_H=GPS_H-interp_file(Lat,Long,Height_filenames.file_AG2020) # AUSGEOID2020 file
+    return [AHD_H]
+
+def DOV(Lat,Long):   
+    # Convert to GPS    
+    DOV_PM=interp_file(Lat,Long,Height_filenames.file_DOV_PM) # AVWS file
+    # Convert to AHD
+    DOV_PV=interp_file(Lat,Long,Height_filenames.file_DOV_PV) # AUSGEOID2020 file
+    return [DOV_PM,DOV_PV]

geodepy/tests/test_heights_s3.py

@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 08 10:52:07 2019
+
+@author: u84157
+"""
+import unittest
+import numpy as np
+import geodepy.heights_s3
+
+#___________________________________________________________________________#
+## Some test Cases to run, check againts ga.gov.au/ausgeoid
+# Test Coordinates
+Long=148.716
+Lat=-25.716
+Height=0
+# What the output should be
+AVWS_H=-40.79437480
+AVWS_H_STD=0.05786271
+AHD_H=-42.231267
+AHD_H_STD=0.10002191
+DOVPM=-13.71323159
+DOVPV=-2.4971921
+
+class TestHeights(unittest.TestCase):
+    def test_AVWS_H(self):
+        self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7)
+    def test_AVWS_H_STD(self):
+        self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7)
+    def test_DOVPV(self):
+        self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.DOV(Lat,Long)[1]),DOVPV,7)
+    def test_DOVPM(self):
+        self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.DOV(Lat,Long)[0]),DOVPM,7)
+    def test_AHD_H(self):
+        self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7)
+    def test_AHD_H_STD(self):
+        self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7)
+
+
+if __name__ == '__main__':
+    unittest.main()

setup.py

@@ -9,4 +9,4 @@
       author_email='geodesy@ga.gov.au',
       license='Apache License 2.0',
       packages=['geodepy'],
-      install_requires=['numpy'])
+      install_requires=['numpy', 'scipy', 'gdal'])