diff --git a/.travis.yml b/.travis.yml index 71c50ef..72313a3 100644 --- a/.travis.yml +++ b/.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 diff --git a/geodepy/Height_filenames.py b/geodepy/Height_filenames.py new file mode 100644 index 0000000..fb58c9e --- /dev/null +++ b/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' \ No newline at end of file diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py new file mode 100644 index 0000000..274f0ea --- /dev/null +++ b/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] diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py new file mode 100644 index 0000000..5a6a7b5 --- /dev/null +++ b/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() diff --git a/setup.py b/setup.py index 3cd47d4..a9e558e 100644 --- a/setup.py +++ b/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'])