From 46434349da1904e35381c223765e65136f2bb25b Mon Sep 17 00:00:00 2001 From: harry093 Date: Fri, 3 Apr 2020 11:23:53 +1100 Subject: [PATCH] Renamed heights_s3 module to height Added contents of Height_filenames to constants module Added contents of Normal_Correction to height module Refactored tests etc. --- {geodepy => Standalone}/ntv2reader.py | 0 geodepy/Height_filenames.py | 13 -- geodepy/Normal_Correction.py | 92 ----------- geodepy/constants.py | 9 + geodepy/geoid.py | 23 --- geodepy/gnss.py | 230 ++++++++++++++++++++++++++ geodepy/height.py | 167 +++++++++++++++++++ geodepy/heights_s3.py | 90 ---------- geodepy/tests/test_NC.py | 33 ---- geodepy/tests/test_height.py | 59 +++++++ geodepy/tests/test_heights_s3.py | 41 ----- 11 files changed, 465 insertions(+), 292 deletions(-) rename {geodepy => Standalone}/ntv2reader.py (100%) delete mode 100644 geodepy/Height_filenames.py delete mode 100644 geodepy/Normal_Correction.py delete mode 100644 geodepy/geoid.py create mode 100644 geodepy/gnss.py create mode 100644 geodepy/height.py delete mode 100644 geodepy/heights_s3.py delete mode 100644 geodepy/tests/test_NC.py create mode 100644 geodepy/tests/test_height.py delete mode 100644 geodepy/tests/test_heights_s3.py diff --git a/geodepy/ntv2reader.py b/Standalone/ntv2reader.py similarity index 100% rename from geodepy/ntv2reader.py rename to Standalone/ntv2reader.py diff --git a/geodepy/Height_filenames.py b/geodepy/Height_filenames.py deleted file mode 100644 index 5639f34..0000000 --- a/geodepy/Height_filenames.py +++ /dev/null @@ -1,13 +0,0 @@ -# -*- 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' -file_GRAV_BA='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/Bouguer_Grav_RELEASE20191107.tif' diff --git a/geodepy/Normal_Correction.py b/geodepy/Normal_Correction.py deleted file mode 100644 index ecec6b6..0000000 --- a/geodepy/Normal_Correction.py +++ /dev/null @@ -1,92 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Spyder Editor - -This is a temporary script file. -""" -import gdal -import numpy as np -import geodepy.Height_filenames as Height_filenames -from scipy.interpolate import griddata - -def mean_normal_grav(Lat,h): - # GRS 80 constants - a=6378137 - b=6356752.3141 - omega=7292115*(10**-11) - e2=0.00669438002290 - GM=3986005*10**8 - k=0.001931851353 - # GRS80 normal gravity - EllGrav=(10**5)*9.7803267715*(1+k*(np.sin(Lat*np.pi/180)**2))/np.sqrt(1-e2*(np.sin(Lat*np.pi/180)**2)) - FA=-((2*(EllGrav/a)*(1+(a-b)/a + omega**2*a**2*b/GM - 2*(a-b)/a*(np.sin(Lat*np.pi/180)**2))*(h**2)/2-3*(EllGrav/a**2)*(h**3)/3)/h) - mean_normal_g=(EllGrav+FA)*(10**-5) - return mean_normal_g - -def normal_grav(Lat,h): - # GRS 80 constants - a=6378137 - b=6356752.3141 - omega=7292115*(10**-11) - e2=0.00669438002290 - GM=3986005*10**8 - k=0.001931851353 - # GRS80 normal gravity - EllGrav=(10**5)*9.7803267715*(1+k*(np.sin(Lat*np.pi/180)**2))/np.sqrt(1-e2*(np.sin(Lat*np.pi/180)**2)) - FA=-(2*EllGrav*h/a)*(1+(a-b)/a+omega**2*a**2*b/GM-2*(a-b)/a*(np.sin(Lat*np.pi/180)**2))+3*(EllGrav*h**2)/(a**2) - normal_g=(EllGrav+FA)*(10**-5) - return normal_g - -def mean_surface_grav(Lat_A,Long_A,H_A,Lat_B,Long_B,H_B): - Surf_Grav_A=interp_grav(Lat_A,Long_A)*(10**-5)+normal_grav(Lat_A,H_A)+0.0419*2.67*H_A*(10**-5) - Surf_Grav_B=interp_grav(Lat_B,Long_B)*(10**-5)+normal_grav(Lat_B,H_B)+0.0419*2.67*H_B*(10**-5) - mean_g=(Surf_Grav_A+Surf_Grav_B)/2 - return mean_g - -def interp_grav(Lat,Long): - # Grid resolution (known) - res=1.0/60 - # open geotiff file - f = gdal.Open(Height_filenames.file_GRAV_BA) - # load band (akin to a variable in dataset) - band = f.GetRasterBand(1) - # get the pixel width, height, etc. - transform = f.GetGeoTransform() - # 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_g=griddata(pos,Surrounding_data_v,(Long,Lat)) - return interp_g - -def normal_correction(Lat_A,Long_A,H_A,Lat_B,Long_B,H_B): - # ellipsoidal gravity at 45 deg. Lat - Gamma_0=9.8061992115 - # Normal Gravity at the point - normal_g_A=mean_normal_grav(Lat_A,H_A) -# print normal_g_A - normal_g_B=mean_normal_grav(Lat_B,H_B) -# print normal_g_B - dn=H_B-H_A - g=mean_surface_grav(Lat_A,Long_A,H_A,Lat_B,Long_B,H_B) - # print g - NC=(dn*(g-Gamma_0)/Gamma_0)+H_A*(normal_g_A-Gamma_0)/Gamma_0-H_B*(normal_g_B-Gamma_0)/Gamma_0 - return NC,g - - - - - diff --git a/geodepy/constants.py b/geodepy/constants.py index e305149..ad1acaa 100644 --- a/geodepy/constants.py +++ b/geodepy/constants.py @@ -371,3 +371,12 @@ def iers2trans(itrf_from, itrf_to, ref_epoch, tx, ty, tz, sc, rx, ry, rz, d_tx, itrf00to88 = iers2trans('ITRF2000', 'ITRF1988', date(1988, 1, 1), 24.7, 11.5, -97.9, 8.95, 0, 0, -0.18, 0.0, -0.6, -1.4, 0.01, 0, 0, 0.02) + +# The locations of files used in the height module +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' +file_GRAV_BA='/vsicurl/https://geoid.s3-ap-southeast-2.amazonaws.com/AVWS/Bouguer_Grav_RELEASE20191107.tif' diff --git a/geodepy/geoid.py b/geodepy/geoid.py deleted file mode 100644 index bd8bb33..0000000 --- a/geodepy/geoid.py +++ /dev/null @@ -1,23 +0,0 @@ -#!/usr/bin/env python3 - -""" -Geoscience Australia - Python Geodesy Package -Geoid Module - -In Development -""" - -import numpy as np -from geodepy.convert import hp2dec - - -# Define test grid points -nvals = np.array([(hp2dec(-31.51), hp2dec(133.48), -8.806), - (hp2dec(-31.51), hp2dec(133.49), -8.743), - (hp2dec(-31.52), hp2dec(133.48), -8.870), - (hp2dec(-31.52), hp2dec(133.49), -8.805)]) - -lat = hp2dec(-31.515996736) -long = hp2dec(133.483540489) - - diff --git a/geodepy/gnss.py b/geodepy/gnss.py new file mode 100644 index 0000000..60640f0 --- /dev/null +++ b/geodepy/gnss.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python3 + +""" +Geoscience Australia - Python Geodesy Package +GNSS Module + +In Development +""" +from numpy import zeros + + +def read_sinex_estimate(file): + """This function reads in the SOLUTION/ESTIMATE block of a SINEX file. It + returns estimate, a list of tuples: + + estimate = [(code, soln, refEpoch, staX, staY, staZ, staX_sd, staY_sd, + staZ_sd[, velX, velY, velZ, velX_sd, velY_sd, velZ_sd])...] + + where: + * code is the stations's 4-character ID + * soln is the segment of the stations's time series + * refEpoch is the epoch of the solution in the form YY:DOY:SSSSS (YY + is the two digit year, DOY is day of year, and SSSSS is the time of + day in seconds + * sta[XYZ] is the station coordinates in the Cartesian reference frame + * sta[XYZ]_sd is the standard deviation of the station coordinates in + the Cartesian reference frame + * vel[XYZ] is the station velocity in the Cartesian reference frame + * vel[XYZ]_sd is the standard deviation of the station velocity in the + Cartesian reference frame + + Velocities are not included in all SINEX files and so are only returned if + present. + + :param file: the input SINEX file + :return: estimate + """ + + # Create data structures and set variables + lines = [] + estimate = [] + velocities = False + go = False + code = '' + soln = '' + epoch = '' + stax = '' + stay = '' + staz = '' + stax_sd = '' + stay_sd = '' + staz_sd = '' + velx = '' + vely = '' + velz = '' + velx_sd = '' + vely_sd = '' + velz_sd = '' + + # Read the SOLUTION/ESTIMATE block into a list and determine if there is + # any velocity information + with open(file) as f: + for line in f: + if line[:18] == '-SOLUTION/ESTIMATE': + break + if go and line[:11] == '*INDEX TYPE': + pass + elif go: + if line[7:10] == 'VEL': + velocities = True + lines.append(line) + if line[:18] == '+SOLUTION/ESTIMATE': + go = True + + for line in lines: + typ = line[7:11] + if typ == 'STAX': + code = line[14:18] + soln = line[23:26].lstrip() + epoch = line[27:39] + stax = float(line[47:68]) + stax_sd = float(line[69:80]) + elif typ == 'STAY': + stay = float(line[47:68]) + stay_sd = float(line[69:80]) + elif typ == 'STAZ': + staz = float(line[47:68]) + staz_sd = float(line[69:80]) + if not velocities: + info = (code, soln, epoch, stax, stay, staz, stax_sd, stay_sd, + staz_sd) + estimate.append(info) + elif typ == 'VELX': + velx = float(line[47:68]) + velx_sd = float(line[69:80]) + elif typ == 'VELY': + vely = float(line[47:68]) + vely_sd = float(line[69:80]) + elif typ == 'VELZ': + velz = float(line[47:68]) + velz_sd = float(line[69:80]) + info = (code, soln, epoch, stax, stay, staz, stax_sd, stay_sd, + staz_sd, velx, vely, velz, velx_sd, vely_sd, velz_sd) + estimate.append(info) + + return estimate + + +def read_sinex_matrix(file): + + """This function reads in the SOLUTION/MATRIX_ESTIMATE block of a SINEX + file. It returns matrix, a list of tuples: + + matrix = [(code, soln, var_x, covar_xy, covar_xz, var_y, covar_yz, + var_z[, var_v_x, covar_v_xy, covar_v_xz, var_v_y, covar_v_yz, + var_v_z])...] + + where: + * code is the stations's 4-character ID + * soln is the segment of the stations's time series + * var_x is the variance in the X coordinate + * covar_xy is the covariance between the X and the Y coordinates + * covar_xz is the covariance between the X and the Z coordinates + * var_y is the variance in the Y coordinate + * covar_yz is the covariance between the Y and the Z coordinates + * var_z is the variance in the Z coordinate + * var_v_x is the variance in the X velocity + * covar_v_xy is the covariance between the X and the Y velocities + * covar_v_xz is the covariance between the X and the Z velocities + * var_v_y is the variance in the Y velocity + * covar_v_yz is the covariance between the Y and the Z velocities + * var_v_z is the variance in the Z velocity + + Velocities are not included in all SINEX files and so their VCV information + is only returned if they are present. + + :param file: the input SINEX file + :return: matrix + """ + + # Read in the codes (station names) and solutions, and check for velocities + data = read_sinex_estimate('GDA2020_RVS.SNX') + code = [] + soln = [] + velocities = False + for station in data: + code.append(station[0]) + soln.append(station[1]) + if len(data[0]) == 15: + velocities = True + + # Read the SOLUTION/MATRIX_ESTIMATE block into a list and determine if the + # matrix is upper or lower triangular + lines = [] + lower_triangular = False + go = False + with open(file) as f: + for line in f: + if line[:25] == '-SOLUTION/MATRIX_ESTIMATE': + break + if go and line[:12] == '*PARA1 PARA2': + pass + elif go: + lines.append(line) + if line[:25] == '+SOLUTION/MATRIX_ESTIMATE': + if line[26] == 'L': + lower_triangular = True + go = True + + # Create an array containing the matrix elements + if velocities: + n = 6 * int(len(code)) + else: + n = 3 * int(len(code)) + element = zeros((n, n)) + matrix = [] + for line in lines: + col = line.rstrip().split() + for i in range(2, len(col)): + element[int(col[0]) - 1][int(col[1]) + i - 3] = float(col[i]) + if velocities: + if lower_triangular: + for i in range(len(code)): + info = (code[i], soln[i], element[6 * i][6 * i], + element[6 * i + 1][6 * i], + element[6 * i + 1][6 * i + 1], + element[6 * i + 2][6 * i], + element[6 * i + 2][6 * i + 1], + element[6 * i + 2][6 * i + 2], + element[6 * i + 3][6 * i + 3], + element[6 * i + 4][6 * i + 3], + element[6 * i + 4][6 * i + 4], + element[6 * i + 5][6 * i + 3], + element[6 * i + 5][6 * i + 4], + element[6 * i + 5][6 * i + 5]) + matrix.append(info) + else: + for i in range(len(code)): + info = (code[i], soln[i], element[6 * i][6 * i], + element[6 * i][6 * i + 1], element[6 * i][6 * i + 2], + element[6 * i + 1][6 * i + 1], + element[6 * i + 1][6 * i + 2], + element[6 * i + 2][6 * i + 2], + element[6 * i + 3][6 * i + 3], + element[6 * i + 3][6 * i + 4], + element[6 * i + 3][6 * i + 5], + element[6 * i + 4][6 * i + 4], + element[6 * i + 4][6 * i + 5], + element[6 * i + 5][6 * i + 5]) + matrix.append(info) + else: + if lower_triangular: + for i in range(len(code)): + info = (code[i], soln[i], element[3 * i][3 * i], + element[3 * i + 1][3 * i], + element[3 * i + 1][3 * i + 1], + element[3 * i + 2][3 * i], + element[3 * i + 2][3 * i + 1], + element[3 * i + 2][3 * i + 2]) + matrix.append(info) + else: + for i in range(len(code)): + info = (code[i], soln[i], element[3 * i][3 * i], + element[3 * i][3 * i + 1], element[3 * i][3 * i + 2], + element[3 * i + 1][3 * i + 1], + element[3 * i + 1][3 * i + 2], + element[3 * i + 2][3 * i + 2]) + matrix.append(info) + + return matrix diff --git a/geodepy/height.py b/geodepy/height.py new file mode 100644 index 0000000..31d21be --- /dev/null +++ b/geodepy/height.py @@ -0,0 +1,167 @@ +#___________________________________________________________________________# +# 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 geodepy.constants as cons +import gdal +import numpy as np +from scipy.interpolate import griddata +#___________________________________________________________________________# +# 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, cons.file_AVWS) # AVWS file + zeta_std=interp_file(Lat, Long, cons.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, cons.file_AVWS) # AVWS file + zeta_std=interp_file(Lat, Long, cons.file_AVWS_STD) # AVWS STD file + 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_file(Lat, Long, cons.file_AG2020) # AUSGEOID2020 file + # Convert to AVWS + Normal_H=GPS_H-interp_file(Lat, Long, cons.file_AVWS) # AVWS file + return [Normal_H] + +def GPS_to_AHD(Lat,Long,GPS_H): + N=interp_file(Lat, Long, cons.file_AG2020) # AUSGEOID2020 file + N_std=interp_file(Lat, Long, cons.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, cons.file_AG2020) # AUSGEOID2020 file + N_std=interp_file(Lat, Long, cons.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, cons.file_AVWS) # AVWS file + # Convert to AHD + AHD_H=GPS_H-interp_file(Lat, Long, cons.file_AG2020) # AUSGEOID2020 file + return [AHD_H] + +def DOV(Lat,Long): + # Convert to GPS + DOV_PM=interp_file(Lat, Long, cons.file_DOV_PM) # AVWS file + # Convert to AHD + DOV_PV=interp_file(Lat, Long, cons.file_DOV_PV) # AUSGEOID2020 file + return [DOV_PM,DOV_PV] + +def mean_normal_grav(Lat,h): + # GRS 80 constants + a=6378137 + b=6356752.3141 + omega=7292115*(10**-11) + e2=0.00669438002290 + GM=3986005*10**8 + k=0.001931851353 + # GRS80 normal gravity + EllGrav=(10**5)*9.7803267715*(1+k*(np.sin(Lat*np.pi/180)**2))/np.sqrt(1-e2*(np.sin(Lat*np.pi/180)**2)) + FA=-((2*(EllGrav/a)*(1+(a-b)/a + omega**2*a**2*b/GM - 2*(a-b)/a*(np.sin(Lat*np.pi/180)**2))*(h**2)/2-3*(EllGrav/a**2)*(h**3)/3)/h) + mean_normal_g=(EllGrav+FA)*(10**-5) + return mean_normal_g + +def normal_grav(Lat,h): + # GRS 80 constants + a=6378137 + b=6356752.3141 + omega=7292115*(10**-11) + e2=0.00669438002290 + GM=3986005*10**8 + k=0.001931851353 + # GRS80 normal gravity + EllGrav=(10**5)*9.7803267715*(1+k*(np.sin(Lat*np.pi/180)**2))/np.sqrt(1-e2*(np.sin(Lat*np.pi/180)**2)) + FA=-(2*EllGrav*h/a)*(1+(a-b)/a+omega**2*a**2*b/GM-2*(a-b)/a*(np.sin(Lat*np.pi/180)**2))+3*(EllGrav*h**2)/(a**2) + normal_g=(EllGrav+FA)*(10**-5) + return normal_g + +def mean_surface_grav(Lat_A,Long_A,H_A,Lat_B,Long_B,H_B): + Surf_Grav_A=interp_grav(Lat_A,Long_A)*(10**-5)+normal_grav(Lat_A,H_A)+0.0419*2.67*H_A*(10**-5) + Surf_Grav_B=interp_grav(Lat_B,Long_B)*(10**-5)+normal_grav(Lat_B,H_B)+0.0419*2.67*H_B*(10**-5) + mean_g=(Surf_Grav_A+Surf_Grav_B)/2 + return mean_g + +def interp_grav(Lat,Long): + # Grid resolution (known) + res=1.0/60 + # open geotiff file + f = gdal.Open(cons.file_GRAV_BA) + # load band (akin to a variable in dataset) + band = f.GetRasterBand(1) + # get the pixel width, height, etc. + transform = f.GetGeoTransform() + # 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_g=griddata(pos,Surrounding_data_v,(Long,Lat)) + return interp_g + +def normal_correction(Lat_A,Long_A,H_A,Lat_B,Long_B,H_B): + # ellipsoidal gravity at 45 deg. Lat + Gamma_0=9.8061992115 + # Normal Gravity at the point + normal_g_A=mean_normal_grav(Lat_A,H_A) +# print normal_g_A + normal_g_B=mean_normal_grav(Lat_B,H_B) +# print normal_g_B + dn=H_B-H_A + g=mean_surface_grav(Lat_A,Long_A,H_A,Lat_B,Long_B,H_B) + # print g + NC=(dn*(g-Gamma_0)/Gamma_0)+H_A*(normal_g_A-Gamma_0)/Gamma_0-H_B*(normal_g_B-Gamma_0)/Gamma_0 + return NC,g diff --git a/geodepy/heights_s3.py b/geodepy/heights_s3.py deleted file mode 100644 index f31e26e..0000000 --- a/geodepy/heights_s3.py +++ /dev/null @@ -1,90 +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" -#___________________________________________________________________________# -# 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): - # 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_NC.py b/geodepy/tests/test_NC.py deleted file mode 100644 index 58f5abd..0000000 --- a/geodepy/tests/test_NC.py +++ /dev/null @@ -1,33 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Dec 06 14:20:43 2019 - -@author: u84157 -""" - - -import unittest -import numpy as np -import geodepy.Normal_Correction - -#___________________________________________________________________________# -## Some test Cases to run, check againts ga.gov.au/ausgeoid -# Test Coordinates -Long1=141.478560 -Lat1=-31.599542 -Height1=368.40 -Long2=141.478560 -Lat2=-31.599542 -Height2=3690.40 -# What the output should be -RECOVERED_GRAV=9.79062606 -NC=0.80179403 - -class TestNC(unittest.TestCase): - def test_NC(self): - self.assertAlmostEqual(np.asscalar(geodepy.Normal_Correction.normal_correction(Lat1,Long1,Height1,Lat2,Long2,Height2)[0]),NC,7) - def test_Grav(self): - self.assertAlmostEqual(np.asscalar(geodepy.Normal_Correction.normal_correction(Lat1,Long1,Height1,Lat2,Long2,Height2)[1]),RECOVERED_GRAV,7) - -if __name__ == '__main__': - unittest.main() \ No newline at end of file diff --git a/geodepy/tests/test_height.py b/geodepy/tests/test_height.py new file mode 100644 index 0000000..529f06b --- /dev/null +++ b/geodepy/tests/test_height.py @@ -0,0 +1,59 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Nov 08 10:52:07 2019 + +@author: u84157 +""" +import unittest +import numpy as np +import geodepy.height + +#___________________________________________________________________________# +## 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.height.GPS_to_AVWS(Lat, Long, Height)[0]), AVWS_H, 7) + def test_AVWS_H_STD(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.GPS_to_AVWS(Lat, Long, Height)[1]), AVWS_H_STD, 7) + def test_DOVPV(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.DOV(Lat, Long)[1]), DOVPV, 7) + def test_DOVPM(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.DOV(Lat, Long)[0]), DOVPM, 7) + def test_AHD_H(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.GPS_to_AHD(Lat, Long, Height)[0]), AHD_H, 7) + def test_AHD_H_STD(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.GPS_to_AHD(Lat, Long, Height)[1]), AHD_H_STD, 7) + +#___________________________________________________________________________# +## Some test Cases to run, check againts ga.gov.au/ausgeoid +# Test Coordinates +Long1=141.478560 +Lat1=-31.599542 +Height1=368.40 +Long2=141.478560 +Lat2=-31.599542 +Height2=3690.40 +# What the output should be +RECOVERED_GRAV=9.79062606 +NC=0.80179403 + +class TestNC(unittest.TestCase): + def test_NC(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.normal_correction(Lat1,Long1,Height1,Lat2,Long2,Height2)[0]),NC,7) + def test_Grav(self): + self.assertAlmostEqual(np.asscalar(geodepy.height.normal_correction(Lat1,Long1,Height1,Lat2,Long2,Height2)[1]),RECOVERED_GRAV,7) + +if __name__ == '__main__': + unittest.main() diff --git a/geodepy/tests/test_heights_s3.py b/geodepy/tests/test_heights_s3.py deleted file mode 100644 index 5a6a7b5..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 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()