From b5b3849d5c33c3c2e7ccfa7bab51811177d9c8eb Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Fri, 8 Nov 2019 11:41:52 +1100 Subject: [PATCH 01/17] Add files via upload --- geodepy/tests/test_heights.py | 41 +++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 geodepy/tests/test_heights.py diff --git a/geodepy/tests/test_heights.py b/geodepy/tests/test_heights.py new file mode 100644 index 0000000..54a9c1f --- /dev/null +++ b/geodepy/tests/test_heights.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 + +#___________________________________________________________________________# +## 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(heights.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7) + def test_AVWS_H_STD(self): + self.assertAlmostEqual(np.asscalar(heights.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7) + def test_DOVPV(self): + self.assertAlmostEqual(np.asscalar(heights.interp_DOVPV(Lat,Long)[0]),DOVPV,7) + def test_DOVPM(self): + self.assertAlmostEqual(np.asscalar(heights.interp_DOVPM(Lat,Long)[0]),DOVPM,7) + def test_AHD_H(self): + self.assertAlmostEqual(np.asscalar(heights.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7) + def test_AHD_H_STD(self): + self.assertAlmostEqual(np.asscalar(heights.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file From 5d0aeeccba3f6a7bcc097bcdeafbcde4b7839449 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Fri, 8 Nov 2019 11:42:24 +1100 Subject: [PATCH 02/17] Add files via upload --- geodepy/heights.py | 233 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 233 insertions(+) create mode 100644 geodepy/heights.py diff --git a/geodepy/heights.py b/geodepy/heights.py new file mode 100644 index 0000000..470649a --- /dev/null +++ b/geodepy/heights.py @@ -0,0 +1,233 @@ +#___________________________________________________________________________# +# 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" +# File names are hard coded +# Some exmaples of how the functions can be used are give towards the end +#___________________________________________________________________________# +# Import dependencies +import gdal +import numpy as np +from scipy.interpolate import griddata +#___________________________________________________________________________# +# Interpolation functions +def interp_DOVPV(Lat,Long): + # Import the DOVPM file + f = gdal.Open('DOV_PV.tif') + # 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_DOVPV_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_DOVPV_val + +def interp_DOVPM(Lat,Long): + # Import the DOVPM file + f = gdal.Open('DOV_PM.tif') + # 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_DOVPM_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_DOVPM_val + +def interp_AUSGeoid2020(Lat,Long): + # open geotiff file + f = gdal.Open('AUSGeoid2020_RELEASEV20170908.tif') + # 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_AUSGeoid2020_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AUSGeoid2020_val + +def interp_AUSGeoid2020_STD(Lat,Long): + # open geotiff file + f = gdal.Open('AUSGeoid2020_RELEASEV20170908_err.tif') + # 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_AUSGeoid2020_STD_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AUSGeoid2020_STD_val + +def interp_AVWS(Lat,Long): + # open geotiff file + f = gdal.Open('AVWS_20191107.tif') + # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AVWS_val + +def interp_AVWS_STD(Lat,Long): + # open geotiff file + f = gdal.Open('AVWS_STD_20191107.tif') + # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AVWS_val + +#___________________________________________________________________________# +# Functions to handle the conversions from one height to another +def GPS_to_AVWS(Lat,Long,GPS_H): + zeta=interp_AVWS(Lat,Long) + zeta_std=interp_AVWS_STD(Lat,Long) + NORMAL_H=GPS_H-zeta + return [NORMAL_H,zeta_std] + +def AVWS_to_GPS(Lat,Long,AVWS_H): + zeta=interp_AVWS(Lat,Long) + zeta_std=interp_AVWS_STD(Lat,Long) + GPS_H=AVWS_H+zeta + return [GPS_H,zeta_std] + +def AHD_to_AVWS(Lat,Long,AHD_H): + # Convert to GPS + GPS_H=AHD_H+interp_AUSGeoid2020(Lat,Long) + # Convert to AVWS + Normal_H=GPS_H-interp_AVWS(Lat,Long) + return [Normal_H] + +def GPS_to_AHD(Lat,Long,GPS_H): + N=interp_AUSGeoid2020(Lat,Long) + N_std=interp_AUSGeoid2020_STD(Lat,Long) + AHD_H=GPS_H-N + return [AHD_H,N_std] + +def AHD_to_GPS(Lat,Long,AHD_H): + N=interp_AUSGeoid2020(Lat,Long) + N_std=interp_AUSGeoid2020_STD(Lat,Long) + 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_AVWS(Lat,Long) + # Convert to AHD + AHD_H=GPS_H-interp_AUSGeoid2020(Lat,Long) + return [AHD_H] + + + + From ca52ad3d0867ad77d8a94832717b1f776dfdde39 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 11 Nov 2019 16:00:02 +1100 Subject: [PATCH 03/17] Delete heights.py --- geodepy/heights.py | 233 --------------------------------------------- 1 file changed, 233 deletions(-) delete mode 100644 geodepy/heights.py diff --git a/geodepy/heights.py b/geodepy/heights.py deleted file mode 100644 index 470649a..0000000 --- a/geodepy/heights.py +++ /dev/null @@ -1,233 +0,0 @@ -#___________________________________________________________________________# -# 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" -# File names are hard coded -# Some exmaples of how the functions can be used are give towards the end -#___________________________________________________________________________# -# Import dependencies -import gdal -import numpy as np -from scipy.interpolate import griddata -#___________________________________________________________________________# -# Interpolation functions -def interp_DOVPV(Lat,Long): - # Import the DOVPM file - f = gdal.Open('DOV_PV.tif') - # 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_DOVPV_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_DOVPV_val - -def interp_DOVPM(Lat,Long): - # Import the DOVPM file - f = gdal.Open('DOV_PM.tif') - # 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_DOVPM_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_DOVPM_val - -def interp_AUSGeoid2020(Lat,Long): - # open geotiff file - f = gdal.Open('AUSGeoid2020_RELEASEV20170908.tif') - # 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_AUSGeoid2020_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AUSGeoid2020_val - -def interp_AUSGeoid2020_STD(Lat,Long): - # open geotiff file - f = gdal.Open('AUSGeoid2020_RELEASEV20170908_err.tif') - # 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_AUSGeoid2020_STD_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AUSGeoid2020_STD_val - -def interp_AVWS(Lat,Long): - # open geotiff file - f = gdal.Open('AVWS_20191107.tif') - # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AVWS_val - -def interp_AVWS_STD(Lat,Long): - # open geotiff file - f = gdal.Open('AVWS_STD_20191107.tif') - # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AVWS_val - -#___________________________________________________________________________# -# Functions to handle the conversions from one height to another -def GPS_to_AVWS(Lat,Long,GPS_H): - zeta=interp_AVWS(Lat,Long) - zeta_std=interp_AVWS_STD(Lat,Long) - NORMAL_H=GPS_H-zeta - return [NORMAL_H,zeta_std] - -def AVWS_to_GPS(Lat,Long,AVWS_H): - zeta=interp_AVWS(Lat,Long) - zeta_std=interp_AVWS_STD(Lat,Long) - GPS_H=AVWS_H+zeta - return [GPS_H,zeta_std] - -def AHD_to_AVWS(Lat,Long,AHD_H): - # Convert to GPS - GPS_H=AHD_H+interp_AUSGeoid2020(Lat,Long) - # Convert to AVWS - Normal_H=GPS_H-interp_AVWS(Lat,Long) - return [Normal_H] - -def GPS_to_AHD(Lat,Long,GPS_H): - N=interp_AUSGeoid2020(Lat,Long) - N_std=interp_AUSGeoid2020_STD(Lat,Long) - AHD_H=GPS_H-N - return [AHD_H,N_std] - -def AHD_to_GPS(Lat,Long,AHD_H): - N=interp_AUSGeoid2020(Lat,Long) - N_std=interp_AUSGeoid2020_STD(Lat,Long) - 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_AVWS(Lat,Long) - # Convert to AHD - AHD_H=GPS_H-interp_AUSGeoid2020(Lat,Long) - return [AHD_H] - - - - From d910c1113cbb83507e63cb297cabdf24d9a0d1bf Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 11 Nov 2019 16:00:23 +1100 Subject: [PATCH 04/17] Add files via upload --- geodepy/heights_s3.py | 234 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) create mode 100644 geodepy/heights_s3.py diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py new file mode 100644 index 0000000..d2e6e95 --- /dev/null +++ b/geodepy/heights_s3.py @@ -0,0 +1,234 @@ +#___________________________________________________________________________# +# 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" +# File names are hard coded +# Some exmaples of how the functions can be used are give towards the end +#___________________________________________________________________________# +# Import dependencies +import gdal +import numpy as np +from scipy.interpolate import griddata +#___________________________________________________________________________# +# Interpolation functions +def interp_DOVPV(Lat,Long): + # Import the DOVPM file + f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/DOV_PV.tif') + # 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_DOVPV_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_DOVPV_val + +def interp_DOVPM(Lat,Long): + # Import the DOVPM file + f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/DOV_PM.tif') + # 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_DOVPM_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_DOVPM_val + +def interp_AUSGeoid2020(Lat,Long): + # open geotiff file + f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AUSGeoid2020_RELEASEV20170908.tif') + # 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_AUSGeoid2020_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AUSGeoid2020_val + +def interp_AUSGeoid2020_STD(Lat,Long): + # open geotiff file + f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AUSGeoid2020_RELEASEV20170908_err.tif') + # 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_AUSGeoid2020_STD_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AUSGeoid2020_STD_val + +def interp_AVWS(Lat,Long): + # open geotiff file + f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AVWS_20191107.tif') + # f = gdal.Open('AVWS_20191107.tif') + # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AVWS_val + +def interp_AVWS_STD(Lat,Long): + # open geotiff file + f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AVWS_STD_20191107.tif') + # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') + return interp_AVWS_val + +#___________________________________________________________________________# +# Functions to handle the conversions from one height to another +def GPS_to_AVWS(Lat,Long,GPS_H): + zeta=interp_AVWS(Lat,Long) + zeta_std=interp_AVWS_STD(Lat,Long) + NORMAL_H=GPS_H-zeta + return [NORMAL_H,zeta_std] + +def AVWS_to_GPS(Lat,Long,AVWS_H): + zeta=interp_AVWS(Lat,Long) + zeta_std=interp_AVWS_STD(Lat,Long) + GPS_H=AVWS_H+zeta + return [GPS_H,zeta_std] + +def AHD_to_AVWS(Lat,Long,AHD_H): + # Convert to GPS + GPS_H=AHD_H+interp_AUSGeoid2020(Lat,Long) + # Convert to AVWS + Normal_H=GPS_H-interp_AVWS(Lat,Long) + return [Normal_H] + +def GPS_to_AHD(Lat,Long,GPS_H): + N=interp_AUSGeoid2020(Lat,Long) + N_std=interp_AUSGeoid2020_STD(Lat,Long) + AHD_H=GPS_H-N + return [AHD_H,N_std] + +def AHD_to_GPS(Lat,Long,AHD_H): + N=interp_AUSGeoid2020(Lat,Long) + N_std=interp_AUSGeoid2020_STD(Lat,Long) + 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_AVWS(Lat,Long) + # Convert to AHD + AHD_H=GPS_H-interp_AUSGeoid2020(Lat,Long) + return [AHD_H] + + + + From 89d205d86399cc0bc6660c93dccacdff62fc963a Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 11 Nov 2019 16:00:53 +1100 Subject: [PATCH 05/17] Delete test_heights.py --- geodepy/tests/test_heights.py | 41 ----------------------------------- 1 file changed, 41 deletions(-) delete mode 100644 geodepy/tests/test_heights.py diff --git a/geodepy/tests/test_heights.py b/geodepy/tests/test_heights.py deleted file mode 100644 index 54a9c1f..0000000 --- a/geodepy/tests/test_heights.py +++ /dev/null @@ -1,41 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Nov 08 10:52:07 2019 - -@author: u84157 -""" -import unittest -import numpy as np -import geodepy.heights - -#___________________________________________________________________________# -## 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(heights.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7) - def test_AVWS_H_STD(self): - self.assertAlmostEqual(np.asscalar(heights.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7) - def test_DOVPV(self): - self.assertAlmostEqual(np.asscalar(heights.interp_DOVPV(Lat,Long)[0]),DOVPV,7) - def test_DOVPM(self): - self.assertAlmostEqual(np.asscalar(heights.interp_DOVPM(Lat,Long)[0]),DOVPM,7) - def test_AHD_H(self): - self.assertAlmostEqual(np.asscalar(heights.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7) - def test_AHD_H_STD(self): - self.assertAlmostEqual(np.asscalar(heights.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) - - -if __name__ == '__main__': - unittest.main() \ No newline at end of file From 11070bffdfc94a6d66a5d45ac1389ed54b252483 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 18 Nov 2019 13:39:54 +1100 Subject: [PATCH 06/17] Add files via upload --- geodepy/tests/test_heights_s3.py | 41 ++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 geodepy/tests/test_heights_s3.py diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py new file mode 100644 index 0000000..2aa16ad --- /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 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(heights_s3.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7) + def test_AVWS_H_STD(self): + self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7) + def test_DOVPV(self): + self.assertAlmostEqual(np.asscalar(heights_s3.interp_DOVPV(Lat,Long)[0]),DOVPV,7) + def test_DOVPM(self): + self.assertAlmostEqual(np.asscalar(heights_s3.interp_DOVPM(Lat,Long)[0]),DOVPM,7) + def test_AHD_H(self): + self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7) + def test_AHD_H_STD(self): + self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file From 21b9de400070b29106a7ac2030a6786296d7075b Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Fri, 29 Nov 2019 09:20:31 +1100 Subject: [PATCH 07/17] Delete heights_s3.py --- geodepy/heights_s3.py | 234 ------------------------------------------ 1 file changed, 234 deletions(-) delete mode 100644 geodepy/heights_s3.py diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py deleted file mode 100644 index d2e6e95..0000000 --- a/geodepy/heights_s3.py +++ /dev/null @@ -1,234 +0,0 @@ -#___________________________________________________________________________# -# 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" -# File names are hard coded -# Some exmaples of how the functions can be used are give towards the end -#___________________________________________________________________________# -# Import dependencies -import gdal -import numpy as np -from scipy.interpolate import griddata -#___________________________________________________________________________# -# Interpolation functions -def interp_DOVPV(Lat,Long): - # Import the DOVPM file - f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/DOV_PV.tif') - # 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_DOVPV_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_DOVPV_val - -def interp_DOVPM(Lat,Long): - # Import the DOVPM file - f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/DOV_PM.tif') - # 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_DOVPM_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_DOVPM_val - -def interp_AUSGeoid2020(Lat,Long): - # open geotiff file - f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AUSGeoid2020_RELEASEV20170908.tif') - # 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_AUSGeoid2020_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AUSGeoid2020_val - -def interp_AUSGeoid2020_STD(Lat,Long): - # open geotiff file - f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AUSGeoid2020_RELEASEV20170908_err.tif') - # 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_AUSGeoid2020_STD_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AUSGeoid2020_STD_val - -def interp_AVWS(Lat,Long): - # open geotiff file - f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AVWS_20191107.tif') - # f = gdal.Open('AVWS_20191107.tif') - # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AVWS_val - -def interp_AVWS_STD(Lat,Long): - # open geotiff file - f = gdal.Open('/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/AVWS_STD_20191107.tif') - # 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_AVWS_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic') - return interp_AVWS_val - -#___________________________________________________________________________# -# Functions to handle the conversions from one height to another -def GPS_to_AVWS(Lat,Long,GPS_H): - zeta=interp_AVWS(Lat,Long) - zeta_std=interp_AVWS_STD(Lat,Long) - NORMAL_H=GPS_H-zeta - return [NORMAL_H,zeta_std] - -def AVWS_to_GPS(Lat,Long,AVWS_H): - zeta=interp_AVWS(Lat,Long) - zeta_std=interp_AVWS_STD(Lat,Long) - GPS_H=AVWS_H+zeta - return [GPS_H,zeta_std] - -def AHD_to_AVWS(Lat,Long,AHD_H): - # Convert to GPS - GPS_H=AHD_H+interp_AUSGeoid2020(Lat,Long) - # Convert to AVWS - Normal_H=GPS_H-interp_AVWS(Lat,Long) - return [Normal_H] - -def GPS_to_AHD(Lat,Long,GPS_H): - N=interp_AUSGeoid2020(Lat,Long) - N_std=interp_AUSGeoid2020_STD(Lat,Long) - AHD_H=GPS_H-N - return [AHD_H,N_std] - -def AHD_to_GPS(Lat,Long,AHD_H): - N=interp_AUSGeoid2020(Lat,Long) - N_std=interp_AUSGeoid2020_STD(Lat,Long) - 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_AVWS(Lat,Long) - # Convert to AHD - AHD_H=GPS_H-interp_AUSGeoid2020(Lat,Long) - return [AHD_H] - - - - From 3d63e930546ff11e4bf32c60812bc8c76051f24f Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Fri, 29 Nov 2019 09:21:30 +1100 Subject: [PATCH 08/17] Add files via upload --- geodepy/Height_filenames.py | 12 +++++ geodepy/heights_s3.py | 92 +++++++++++++++++++++++++++++++++++++ 2 files changed, 104 insertions(+) create mode 100644 geodepy/Height_filenames.py create mode 100644 geodepy/heights_s3.py 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..04249be --- /dev/null +++ b/geodepy/heights_s3.py @@ -0,0 +1,92 @@ +#___________________________________________________________________________# +# 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" +# File names are hard coded +# Some exmaples of how the functions can be used are give towards the end +#___________________________________________________________________________# +# Import dependencies +import gdal +import numpy as np +from scipy.interpolate import griddata +import 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] From 190df2f730455d2cb0cd390b0fd743995ec4f624 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Fri, 29 Nov 2019 09:21:52 +1100 Subject: [PATCH 09/17] Delete test_heights_s3.py --- geodepy/tests/test_heights_s3.py | 41 -------------------------------- 1 file changed, 41 deletions(-) delete mode 100644 geodepy/tests/test_heights_s3.py diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py deleted file mode 100644 index 2aa16ad..0000000 --- a/geodepy/tests/test_heights_s3.py +++ /dev/null @@ -1,41 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Nov 08 10:52:07 2019 - -@author: u84157 -""" -import unittest -import numpy as np -import 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(heights_s3.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7) - def test_AVWS_H_STD(self): - self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7) - def test_DOVPV(self): - self.assertAlmostEqual(np.asscalar(heights_s3.interp_DOVPV(Lat,Long)[0]),DOVPV,7) - def test_DOVPM(self): - self.assertAlmostEqual(np.asscalar(heights_s3.interp_DOVPM(Lat,Long)[0]),DOVPM,7) - def test_AHD_H(self): - self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7) - def test_AHD_H_STD(self): - self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) - - -if __name__ == '__main__': - unittest.main() \ No newline at end of file From a1906e7cc6d2d73459c959fe420edf3d7094d89a Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Fri, 29 Nov 2019 09:23:23 +1100 Subject: [PATCH 10/17] Add files via upload --- geodepy/tests/test_heights_s3.py | 41 ++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 geodepy/tests/test_heights_s3.py diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py new file mode 100644 index 0000000..4fce0ca --- /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 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(heights_s3.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7) + def test_AVWS_H_STD(self): + self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7) + def test_DOVPV(self): + self.assertAlmostEqual(np.asscalar(heights_s3.DOV(Lat,Long)[1]),DOVPV,7) + def test_DOVPM(self): + self.assertAlmostEqual(np.asscalar(heights_s3.DOV(Lat,Long)[0]),DOVPM,7) + def test_AHD_H(self): + self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7) + def test_AHD_H_STD(self): + self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file From 1b1d748ad57410b8a76ebb9bae6e7e794bb59a3a Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 15:22:24 +1100 Subject: [PATCH 11/17] Update geodepy/tests/test_heights_s3.py Co-Authored-By: Josh Batchelor <35472516+BatchelorJ@users.noreply.github.com> --- geodepy/tests/test_heights_s3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py index 4fce0ca..8d0b181 100644 --- a/geodepy/tests/test_heights_s3.py +++ b/geodepy/tests/test_heights_s3.py @@ -6,7 +6,7 @@ """ import unittest import numpy as np -import heights_s3 +import geodepy.heights_s3 #___________________________________________________________________________# ## Some test Cases to run, check againts ga.gov.au/ausgeoid @@ -38,4 +38,4 @@ def test_AHD_H_STD(self): if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() From a8c7186304aa237b040b0766c3e937ac2933ae80 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 15:25:50 +1100 Subject: [PATCH 12/17] Update heights_s3.py --- geodepy/heights_s3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py index 04249be..6997734 100644 --- a/geodepy/heights_s3.py +++ b/geodepy/heights_s3.py @@ -12,7 +12,7 @@ import gdal import numpy as np from scipy.interpolate import griddata -import Height_filenames +import geodepy.Height_filenames #___________________________________________________________________________# # Interpolation functions def interp_file(Lat,Long,file): From 1924737947925d6c58e8eff785ed3c72efca33ad Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 15:37:38 +1100 Subject: [PATCH 13/17] Update heights_s3.py --- geodepy/heights_s3.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py index 6997734..9ec5700 100644 --- a/geodepy/heights_s3.py +++ b/geodepy/heights_s3.py @@ -5,8 +5,6 @@ # AWVS/AHD and Vice Versa # Gridded data used for the varisous reference surfaces are geotiff files # These allow direct access remotely using "gdal" -# File names are hard coded -# Some exmaples of how the functions can be used are give towards the end #___________________________________________________________________________# # Import dependencies import gdal From ef09ae20f8e61083f58b84854214f33f47336ab3 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 16:04:52 +1100 Subject: [PATCH 14/17] Add gdal and scipy dependencies --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From caa0c23fb2c3a9f4fc70c0e1df935a39b64675ee Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 16:08:49 +1100 Subject: [PATCH 15/17] Update test_heights_s3.py --- geodepy/tests/test_heights_s3.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py index 8d0b181..5a6a7b5 100644 --- a/geodepy/tests/test_heights_s3.py +++ b/geodepy/tests/test_heights_s3.py @@ -24,17 +24,17 @@ class TestHeights(unittest.TestCase): def test_AVWS_H(self): - self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AVWS(Lat,Long,Height)[0]),AVWS_H,7) + 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(heights_s3.GPS_to_AVWS(Lat,Long,Height)[1]),AVWS_H_STD,7) + 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(heights_s3.DOV(Lat,Long)[1]),DOVPV,7) + self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.DOV(Lat,Long)[1]),DOVPV,7) def test_DOVPM(self): - self.assertAlmostEqual(np.asscalar(heights_s3.DOV(Lat,Long)[0]),DOVPM,7) + self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.DOV(Lat,Long)[0]),DOVPM,7) def test_AHD_H(self): - self.assertAlmostEqual(np.asscalar(heights_s3.GPS_to_AHD(Lat,Long,Height)[0]),AHD_H,7) + 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(heights_s3.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) + self.assertAlmostEqual(np.asscalar(geodepy.heights_s3.GPS_to_AHD(Lat,Long,Height)[1]),AHD_H_STD,7) if __name__ == '__main__': From 4b41c91bd0f4b78d7e0a5f35299fbf3e2aef4628 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 16:13:44 +1100 Subject: [PATCH 16/17] Update heights_s3.py --- geodepy/heights_s3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py index 9ec5700..274f0ea 100644 --- a/geodepy/heights_s3.py +++ b/geodepy/heights_s3.py @@ -10,7 +10,7 @@ import gdal import numpy as np from scipy.interpolate import griddata -import geodepy.Height_filenames +import geodepy.Height_filenames as Height_filenames #___________________________________________________________________________# # Interpolation functions def interp_file(Lat,Long,file): From 3c9a2f4c46f368d040b1fc65a69af3c0dc803bd2 Mon Sep 17 00:00:00 2001 From: JackMcCubbineGA <48938190+JackMcCubbineGA@users.noreply.github.com> Date: Mon, 2 Dec 2019 16:20:56 +1100 Subject: [PATCH 17/17] Update setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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'])