From 29a0aa3c12f83ed3bbc84ea73bc5a0c30250dabe Mon Sep 17 00:00:00 2001 From: Jaimie Dodd Date: Mon, 17 Dec 2018 09:58:44 +1000 Subject: [PATCH] NTv2 .gsb file reader & interpolator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Python script reads the binary NTv2 file for AUSGeoid2020, finds the four surrounding nodes for the coordinates supplied and performs bilinear or bicubic interpolation (based on user input, or bicubic if not specified) . The output is a tuple of the ellipsoid to AHD separation, deflection in the prime meridian and deflection in the prime vertical values in that order (NAHD, ξ, η). The script contains 2 classes with multiple functions, but only NTv2ReaderBinary().ntv2reader(file, lat, lon, interpolation_method) needs to be called, where: • file is the complete file location as a string • lat is the latitude of the point of interest in decimal degrees • lon is the longitude of the point of interest in decimal degrees • interpolation_method is either ‘bilinear’ or ‘bicubic’ (by default it is set to ‘bicubic’, so this parameter isn’t mandatory). The script also contains 2 stand alone functions - one for bilinear and one for bicubic interpolation. Signed-off-by: Jaimie Dodd --- geodepy/ntv2reader.py | 621 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 621 insertions(+) create mode 100644 geodepy/ntv2reader.py diff --git a/geodepy/ntv2reader.py b/geodepy/ntv2reader.py new file mode 100644 index 0000000..ce7bcbe --- /dev/null +++ b/geodepy/ntv2reader.py @@ -0,0 +1,621 @@ +import struct +from datetime import datetime as dt +import numpy as np +import pandas as pd + +''' +Created 12/11/2018 + +@Author: Jaimie Dodd + +This script defines classes for reading the AUSGeoid2020 binary NTv2 file, and interpolating the Ellipsoid to AHD +Separation, Deflection in the Prime Meridian and Deflection in the Prime Vertical of a given latitude and longitude via +bilinear or bicubic interpolation (user specified input or bicubic by default). + +Values have been tested against the output of DynAdjust. +''' + +class NTv2ReaderBinary(object): + def __init__(self): + self.numOrec = 0 + self.numSrec = 0 + self.numFile = 0 + self.gsType = '' + self.version = '' + self.systemF = '' + self.systemT = '' + self.majorF = 0 + self.minorF = 0 + self.majorT = 0 + self.minorT = 0 + self.subgrids = pd.DataFrame(columns=['CONTAINS', 'SUB_NAME', 'PARENT', 'CREATED', 'UPDATED', 'S_LAT', 'N_LAT', + 'E_LONG', 'W_LONG', 'LAT_INC', 'LONG_INC', 'GS_COUNT', + 'ELLIPSOID_TO_AHD_SEPARATION', 'DEFLECTION_PRIME_MERIDIAN', + 'DEFLECTION_PRIME_VERTICAL']) + self.ellps2ahdsep = 0 + self.deflectionprimemeridian = 0 + self.deflectionprimevertical = 0 + + def ntv2reader(self, file, lat, lon, interpolation_method='bicubic'): + """ + Function for reading in the binary NTv2 file and + + :param file: NTv2 file (full location as string) + :param lat: latitude of point of interest (decimal degrees) + :param lon: longitude of point of interest (decimal degrees) + :param interpolation_method: 'bilinear' or 'bicubic' ('bicubic' by default) + :return: tuple containing ellipsoid to AHD separation, deflection in prime meridian and deflection in prime + vertical. + + """ + + # If interpolation method specified is something other than bilinear or bicubic, then raise ValueError + interpolation_options = ['bilinear', 'bicubic'] + if interpolation_method not in interpolation_options: + raise ValueError("Invalid Interpolation method. Expected one of: %s" % interpolation_options) + + # Read the NTv2 file as binary + f = open(file, 'rb') + + # NUM_OREC + f.seek(8, 1) + + byte = f.read(4) + self.numOrec = int.from_bytes(byte, byteorder='little') + + f.seek(4, 1) + + # NUM_SREC + f.seek(8, 1) + + byte = f.read(4) + self.numSrec = int.from_bytes(byte, byteorder='little') + + f.seek(4, 1) + + # NUM_FILE + f.seek(8, 1) + + byte = f.read(4) + self.numFile = int.from_bytes(byte, byteorder='little') + + f.seek(4, 1) + + # GS_TYPE + f.seek(8, 1) + + byte = f.read(8) + self.gsType = byte.decode('utf8').strip('\x00').strip() + + # VERSION + f.seek(8, 1) + + byte = f.read(8) + self.version = byte.decode('utf8').strip('\x00').strip() + + # SYSTEM_F + f.seek(8, 1) + + byte = f.read(8) + self.systemF = byte.decode('utf8').strip('\x00').strip() + + # SYSTEM_T + f.seek(8, 1) + + byte = f.read(8) + self.systemT = byte.decode('utf8').strip('\x00').strip() + + # MAJOR_F + f.seek(8, 1) + + byte = f.read(8) + self.majorF = struct.unpack('d', byte)[0] + + # MINOR_F + f.seek(8, 1) + + byte = f.read(8) + self.minorF = struct.unpack('d', byte)[0] + + # MAJOR_T + f.seek(8, 1) + + byte = f.read(8) + self.majorT = struct.unpack('d', byte)[0] + + # MINOR_T + f.seek(8, 1) + + byte = f.read(8) + self.minorT = struct.unpack('d', byte)[0] + + # Convert lat and lon to seconds + if self.gsType == 'SECONDS': + lat = lat * 3600 + lon = lon * -3600 + + # Sub Grids + for i in range(0, self.numFile): + self.subgrids = self.subgrids.append(SubGrid().findsubgrid(f, lat, lon, interpolation_method), + ignore_index=True) + + # Close the NTv2 file + f.close() + + # Filter subgrids dataframe so that only subgrids with values remain + self.subgrids = self.subgrids[self.subgrids['CONTAINS']] + + # If more than one subgrid which contains coordinates, then filter dataframe so subgrid with smallest gridwidth + # remains + if len(self.subgrids) > 1: + self.subgrids = self.subgrids.loc[[self.subgrids.loc[:, 'LAT_INC'].idxmin()]] + + # If subgrids dataframe is not empty, then return values for ellipsoid to AHD separation, deflection of prime + # meridian and deflection of prime vertical + if not self.subgrids.empty: + self.ellps2ahdsep = self.subgrids.iloc[0, -3] + self.deflectionprimemeridian = self.subgrids.iloc[0, -2] + self.deflectionprimevertical = self.subgrids.iloc[0, -1] + return self.ellps2ahdsep, self.deflectionprimemeridian, self.deflectionprimevertical + # If no subgrids contain coordinates, then raise ValueError + else: + raise ValueError("The coordinates supplied are outside the extents of the grid.") + + +class SubGrid(object): + + def __init__(self): + self.subName = '' + self.parent = '' + self.created = '' + self.updated = '' + self.sLat = 0 + self.nLat = 0 + self.eLong = 0 + self.wLong = 0 + self.latInc = 0 + self.longInc = 0 + self.gsCount = 0 + self.Nval = 0 + self.Xi = 0 + self.Eta = 0 + self.contains = False + + def findsubgrid(self, f, lat, lon, interpolation_method='bicubic'): + """ + Function to pull out metadata of a subgrid in the AUSGeoid2020 NTv2 file + + :param f: variable for open file NTv2 being read as binary + :param lat: latitude of the point of interest (seconds) + :param lon: longitude of the point of interest (negative seconds) + :param interpolation_method: interpolation method as specified in ntv2reader() - + 'bilinear' or 'bicubic' ('bicubic' by default) + + :return: subgrid metadata in form of a dictionary or results from bilinear() or bicubic() (depending on + interpolation method specified) if point lies within subgrid. + """ + # SUB_NAME + f.seek(8, 1) + + byte = f.read(8) + self.subName = byte.decode('utf').strip('\x00').strip() + + # PARENT + f.seek(8, 1) + + byte = f.read(8) + self.parent = byte.decode('utf').strip('\x00').strip() + + # CREATED + f.seek(8, 1) + + byte = f.read(8) + self.created = dt.strptime(byte.decode('utf').strip('\x00').strip(), '%d%m%Y').strftime('%d/%m/%Y') + + # UPDATED + f.seek(8, 1) + + byte = f.read(8) + self.updated = dt.strptime(byte.decode('utf').strip('\x00').strip(), '%d%m%Y').strftime('%d/%m/%Y') + + # S_LAT + f.seek(8, 1) + + byte = f.read(8) + self.sLat = round(struct.unpack('d', byte)[0], 3) + + # N_LAT + f.seek(8, 1) + + byte = f.read(8) + self.nLat = round(struct.unpack('d', byte)[0], 3) + + # E_LONG + f.seek(8, 1) + + byte = f.read(8) + self.eLong = round(struct.unpack('d', byte)[0], 3) + + # W_LONG + f.seek(8, 1) + + byte = f.read(8) + self.wLong = round(struct.unpack('d', byte)[0], 3) + + # LAT_INC + f.seek(8, 1) + + byte = f.read(8) + self.latInc = round(struct.unpack('d', byte)[0], 6) + + # LONG_INC + f.seek(8, 1) + + byte = f.read(8) + self.longInc = round(struct.unpack('d', byte)[0], 6) + + # GS_COUNT + f.seek(8, 1) + + byte = f.read(4) + self.gsCount = int.from_bytes(byte, byteorder='little') + + f.seek(4, 1) + + # Check if coordinates fall within subgrid. If yes, return self.contains = True + if self.sLat <= lat < self.nLat and self.eLong <= lon < self.wLong: + self.contains = True + + # If self.contains is True, determine number of columns in grid, and row and column of node to bottom right of + # point of interest, then call relevant interpolation method function + if self.contains is True: + # Determine number of columns + numcols = 1 + int((self.wLong - self.eLong) / self.longInc) + + # Determine row and col numbers of node below right of point + row = int((lat - self.sLat) / self.latInc) + col = int((lon - self.eLong) / self.longInc) + + # If interpolation_method == 'bilinear', call bilinear function + if interpolation_method == 'bilinear': + return self.bilinear(f, lat, lon, numcols, row, col) + # If interpolation_method == 'bicubic', call bicubic function + elif interpolation_method == 'bicubic': + return self.bicubic(f, lat, lon, numcols, row, col) + + # If point is not in subgrid, skip to end of subgrid and return sub grid metadata in form of dictionary + else: + f.seek(16 * self.gsCount, 1) + + return {'CONTAINS': self.contains, 'SUB_NAME': self.subName, 'PARENT': self.parent, 'CREATED': self.created, + 'UPDATED': self.updated, 'S_LAT': self.sLat, 'N_LAT': self.nLat, 'E_LONG': self.eLong, + 'W_LONG': self.wLong, 'LAT_INC': self.latInc, 'LONG_INC': self.longInc, 'GS_COUNT': self.gsCount, + 'ELLIPSOID_TO_AHD_SEPARATION': np.nan, 'DEFLECTION_PRIME_MERIDIAN': np.nan, + 'DEFLECTION_PRIME_VERTICAL': np.nan} + + def bilinear(self, f, lat, lon, numcols, row, col): + + """ + Function to perform bilinear interpolatoin of the Ellipsoid to AHD Separation, Deflection in Prime Meridian and + Deflection in Prime Vertical of a point of interest. + + :param f: variable for open file NTv2 being read as binary + :param lat: latitude of the point of interest (seconds) + :param lon: longitude of the point of interest (negative seconds) + :param numcols: number of columns in grid + :param row: row number of point to the bottom right of the point of interest + :param col: column number of point to the bottom right of the point of interest + + :return: dictionary of subgrid metadata and values for Ellipsoid to AHD Separation, Deflection in Prime Meridian + and Deflection in Prime Vertical of point of interest found via bilinear interpolation. + """ + + # | | + # --o-----o-- + # |4 |3 + # | *P | + # --o-----o-- + # |2 |1 + # + # o - Node position. + + # Determine position in the file of the four surrounding nodes + pos1 = row * numcols + col + pos2 = pos1 + 1 + pos3 = pos1 + numcols + pos4 = pos3 + 1 + + # Navigate to start of posA node + f.seek(16 * pos1, 1) + + # Read in values for nodes A and B + (pos1nval, pos1xi, pos1eta) = self.read_node(f) + (pos2nval, pos2xi, pos2eta) = self.read_node(f) + + # Navigate to beginning of node C + f.seek(16*(pos3 - pos2 - 1), 1) + + # Read in values for nodes C and D + (pos3nval, pos3xi, pos3eta) = self.read_node(f) + (pos4nval, pos4xi, pos4eta) = self.read_node(f) + + # Determine latitude and longitude of node A + lat1 = self.sLat + row * self.latInc + long1 = self.eLong + col * self.longInc + + # Determine interpolation scale factors + x = round((lon - long1) / self.longInc, 6) + y = round((lat - lat1) / self.latInc, 6) + + # Call bilinear interpolation function to determine N value, Xi and Eta + # (Ellipsoid to AHD separation, deflection in prime meridian and deflection in prime vertical). + self.Nval = round(bilinear_interpolation(pos1nval, pos2nval, pos3nval, pos4nval, x, y), 3) + self.Xi = round(bilinear_interpolation(pos1xi, pos2xi, pos3xi, pos4xi, x, y), 2) + self.Eta = round(bilinear_interpolation(pos1eta, pos2eta, pos3eta, pos4eta, x, y), 2) + + # Navigate to the end of the subgrid + f.seek(16 * (self.gsCount - pos4 - 1), 1) + + # Return dictionary of information + return {'CONTAINS': self.contains, 'SUB_NAME': self.subName, 'PARENT': self.parent, 'CREATED': self.created, + 'UPDATED': self.updated, 'S_LAT': self.sLat, 'N_LAT': self.nLat, 'E_LONG': self.eLong, + 'W_LONG': self.wLong, 'LAT_INC': self.latInc, 'LONG_INC': self.longInc, 'GS_COUNT': self.gsCount, + 'ELLIPSOID_TO_AHD_SEPARATION': self.Nval, 'DEFLECTION_PRIME_MERIDIAN': self.Xi, + 'DEFLECTION_PRIME_VERTICAL': self.Eta} + + def bicubic(self, f, lat, lon, numcols, row, col): + """ + Function to perform bicubic interpolatoin of the Ellipsoid to AHD Separation, Deflection in Prime Meridian and + Deflection in Prime Vertical of a point of interest. + + :param f: variable for open file NTv2 being read as binary + :param lat: latitude of the point of interest (seconds) + :param lon: longitude of the point of interest (negative seconds) + :param numcols: number of columns in grid + :param row: row number of point to the bottom right of the point of interest + :param col: column number of point to the bottom right of the point of interest + + :return: dictionary of subgrid metadata and values for Ellipsoid to AHD Separation, Deflection in Prime Meridian + and Deflection in Prime Vertical of point of interest found via bicubic interpolation. + + """ + + # | | | | + # --o-----o-----o-----o-- + # |11 |12 |13 |14 + # | | | | + # --o-----o-----o-----o-- + # |10 |3 |4 |15 + # | | *P | | + # --o-----o-----o-----o-- + # |9 |2 |1 |16 + # | | | | + # --o-----o-----o-----o-- + # |8 |7 |6 |5 + # + # o - Node position. + + # Determine position in the file of the sixteen surrounding nodes + pos1 = row * numcols + col + pos2 = pos1 + 1 + pos3 = pos2 + numcols + pos4 = pos3 - 1 + pos5 = pos4 - 2 * numcols - 1 + pos6 = pos5 + 1 + pos7 = pos6 + 1 + pos8 = pos7 + 1 + pos9 = pos8 + numcols + pos10 = pos9 + numcols + pos11 = pos10 + numcols + pos12 = pos11 - 1 + pos13 = pos12 - 1 + pos14 = pos13 - 1 + pos15 = pos14 - numcols + pos16 = pos15 - numcols + + # Navigate to start of posA node + f.seek(16 * pos5, 1) + + # Read in values for nodes 5-8 + (pos5nval, pos5xi, pos5eta) = self.read_node(f) + (pos6nval, pos6xi, pos6eta) = self.read_node(f) + (pos7nval, pos7xi, pos7eta) = self.read_node(f) + (pos8nval, pos8xi, pos8eta) = self.read_node(f) + + # Navigate to start of pos16 node + f.seek(16 * (pos16 - pos8 - 1), 1) + + # Read in values for nodes 16, 1, 2, and 9 + (pos16nval, pos16xi, pos16eta) = self.read_node(f) + (pos1nval, pos1xi, pos1eta) = self.read_node(f) + (pos2nval, pos2xi, pos2eta) = self.read_node(f) + (pos9nval, pos9xi, pos9eta) = self.read_node(f) + + # Navigate to start of pos15 node + f.seek(16 * (pos15 - pos9 - 1), 1) + + # Read in values for nodes 15, 3, 4 and 10 + (pos15nval, pos15xi, pos15eta) = self.read_node(f) + (pos4nval, pos4xi, pos4eta) = self.read_node(f) + (pos3nval, pos3xi, pos3eta) = self.read_node(f) + (pos10nval, pos10xi, pos10eta) = self.read_node(f) + + # Navigate to start of pos14 node + f.seek(16 * (pos14 - pos10 - 1), 1) + + # Read in values for nodes 11, 12, 13 and 14 + (pos14nval, pos14xi, pos14eta) = self.read_node(f) + (pos13nval, pos13xi, pos13eta) = self.read_node(f) + (pos12nval, pos12xi, pos12eta) = self.read_node(f) + (pos11nval, pos11xi, pos11eta) = self.read_node(f) + + # Determine latitude and longitude of node A + lat1 = self.sLat + row * self.latInc + long1 = self.eLong + col * self.longInc + + # Determine interpolation scale factors + x = round((lon - long1) / self.longInc, 6) + y = round((lat - lat1) / self.latInc, 6) + + # Call bicubic_interpolation to determine nVal, Xi and Eta + # (Ellipsoid to AHD separation, deflection in prime meridian and deflection in prime vertical). + self.Nval = round(bicubic_interpolation(pos1nval, pos2nval, pos3nval, pos4nval, pos5nval, pos6nval, pos7nval, + pos8nval, pos9nval, pos10nval, pos11nval, pos12nval, pos13nval, + pos14nval, pos15nval, pos16nval, x, y), 3) + self.Xi = round(bicubic_interpolation(pos1xi, pos2xi, pos3xi, pos4xi, pos5xi, pos6xi, pos7xi, pos8xi, pos9xi, + pos10xi, pos11xi, pos12xi, pos13xi, pos14xi, pos15xi, pos16xi, x, y), 2) + self.Eta = round(bicubic_interpolation(pos1eta, pos2eta, pos3eta, pos4eta, pos5eta, pos6eta, pos7eta, pos8eta, + pos9eta, pos10eta, pos11eta, pos12eta, pos13eta, pos14eta, pos15eta, + pos16eta, x, y), 2) + + # Navigate to the end of the subgrid + f.seek(16 * (self.gsCount - pos11 - 1), 1) + + # Return dictionary of information + return {'CONTAINS': self.contains, 'SUB_NAME': self.subName, 'PARENT': self.parent, 'CREATED': self.created, + 'UPDATED': self.updated, 'S_LAT': self.sLat, 'N_LAT': self.nLat, 'E_LONG': self.eLong, + 'W_LONG': self.wLong, 'LAT_INC': self.latInc, 'LONG_INC': self.longInc, 'GS_COUNT': self.gsCount, + 'ELLIPSOID_TO_AHD_SEPARATION': self.Nval, 'DEFLECTION_PRIME_MERIDIAN': self.Xi, + 'DEFLECTION_PRIME_VERTICAL': self.Eta} + + # Define function to read in node values + @staticmethod + def read_node(f): + """ + Function to read in values of nodes + + :param f: variable for open file NTv2 being read as binary + :return: Ellipsoid to AHD Separation, Deflection in Prime Meridian and Deflection in Prime Vertical for node. + + """ + # Read ellipsoid to AHD separation value (N) + byte = f.read(4) + nval = round(struct.unpack('f', byte)[0], 3) + # Read deflection in prime meridian value (Xi) + byte = f.read(4) + xi = round(struct.unpack('f', byte)[0], 3) + # Read deflection in prime vertical value (Eta) + byte = f.read(4) + eta = round(struct.unpack('f', byte)[0], 3) + # Skip to beginning of next node + f.seek(4, 1) + # Return values + return nval, xi, eta + + +# Define function for bilinear interpolation +def bilinear_interpolation(n1, n2, n3, n4, x, y): + """ + Bilinear interpolation of value for point of interest (P). + + :param n1: value at node 1 + :param n2: value at node 2 + :param n3: value at node 3 + :param n4: value at node 4 + :param x: interpolation scale factor for x axis + :param y: interpolation scale factor for y axis + + :return: value at node P + + """ + + a0 = n1 + a1 = round(n2 - n1, 3) + a2 = round(n3 - n1, 3) + a3 = round(n1 + n4 - n2 - n3, 3) + p = a0 + (a1 * x) + (a2 * y) + (a3 * x * y) + return p + + +# Define function for performing bicubic interpolation +def bicubic_interpolation(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16, x, y): + """ + Bicubic interpolation of value for point of interest (P). + + :param n1: value at node 1 + :param n2: value at node 2 + :param n3: value at node 3 + :param n4: value at node 4 + :param n5: value at node 5 + :param n6: value at node 6 + :param n7: value at node 7 + :param n8: value at node 8 + :param n9: value at node 9 + :param n10: value at node 10 + :param n11: value at node 11 + :param n12: value at node 12 + :param n13: value at node 13 + :param n14: value at node 14 + :param n15: value at node 16 + :param n16: value at node 17 + :param x: interpolation scale factor for x axis + :param y: interpolation scale factor for y axis + + :return: value at node P + """ + + # Define the inverse of the coefficient matrix + cinv = np.array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [-3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0], + [2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1], + [0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1], + [-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0], + [9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2], + [-6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2], + [2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0], + [-6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1], + [4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1]]) + + # Define x parameters + # Function values + x1 = n1 + x2 = n2 + x3 = n3 + x4 = n4 + # X Derivatives + x5 = round((n2 - n16) / 2, 4) + x6 = round((n9 - n1) / 2, 4) + x7 = round((n10 - n4) / 2, 4) + x8 = round((n3 - n15) / 2, 4) + # Y Derivatives + x9 = round((n4 - n6) / 2, 4) + x10 = round((n3 - n7) / 2, 4) + x11 = round((n12 - n2) / 2, 4) + x12 = round((n13 - n1) / 2, 4) + # Cross Derivatives (XY) + x13 = round((n3 - n7 - n15 + n5) / 4, 4) + x14 = round((n10 - n8 - n4 + n6) / 4, 4) + x15 = round((n11 - n9 - n13 + n1) / 4, 4) + x16 = round((n12 - n2 - n14 + n16) / 4, 4) + + # Create array from x parameters + xarr = np.array([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16]) + + # Multiply the inverse of the coefficient matrix by the array of x values to give array of alpha values + alpha = np.matmul(cinv, xarr) + + # Calculate value at the point of interest + n_p = 0 + for i in range(0, 4): + for j in range(0, 4): + n_p = n_p + alpha[i * 4 + j] * x ** i * y ** j + + # Return the value + return n_p + + +# TEST OF SCRIPT +# Specify AUSGeoid2020 binary NTv2 file location +ntv2file = "C://Git/Python/NTv2.git/AUSGeoid2020_20180201.gsb" +# Assign class to variable +grids = NTv2ReaderBinary() +# Call ntv2reader to determine values at specific point +values = grids.ntv2reader(ntv2file, -26.948643, 145.6548721, 'bilinear') +# Print values +print(values)