Add First Velocity Correction Functionality

geodepy/survey.py

@@ -5,10 +5,50 @@
 Survey Module
 """
 
-from math import sqrt, sin, cos, asin, radians, degrees
+from math import sqrt, sin, cos, asin, radians, degrees, exp
 from geodepy.convert import hp2dec, dec2hp
 
 
+def first_vel_params(wavelength, temp=12, pressure=1013.25, rel_humidity=60):
+    """
+    Calculates the constant First Velocity Correction Parameters C and D
+    for a given set of standard instrument settings
+    :param wavelength: Instrument Carrier Wavelength (micrometres)
+    :param temp: Standard Temperature (degrees Celsius)
+    :param pressure: Standard Pressure (hectopascals or millibars)
+    :param rel_humidity: Standard Relative Humidity (percentage)
+    :return: First Velocity Correction Parameters C and D
+    """
+    group_refractivity = 287.6155 + (4.8866 / wavelength ** 2) + (0.068 / wavelength ** 4)
+    param_d = (273.15 / 1013.25) * group_refractivity
+    saturation_pressure = ((1.0007 + ((3.46 * pressure) * (10 ** -6)))
+                           * 6.1121 * exp((17.502 * temp) / (240.94 + temp)))
+    param_c = ((((group_refractivity * pressure) / (273.15 + temp)) * (273.15 / 1013.25))
+               - (11.27 * ((saturation_pressure * (rel_humidity / 100)) / (273.15 + temp))))
+    return param_c, param_d
+
+
+def first_vel_corrn(dist, first_vel_param, temp, pressure, rel_humidity):
+    """
+    Carries out First Velocity Correction of Electronic Distance Measurement,
+    given correction parameters and atmospheric observations
+    :param dist: Uncorrected Observed Slope Distance
+    :param first_vel_param: Tuple of First Velocity Parameters C and D (see function first_vel_params)
+    :param temp: Observed Temperature (degrees Celsius)
+    :param pressure: Observed Pressure (hectopascals or millibars)
+    :param rel_humidity: Observed Relative Humidity (percentage)
+    :return: Slope Distance with First Velocity Correction applied
+    """
+    param_c = first_vel_param[0]
+    param_d = first_vel_param[1]
+    part_h2o_press = exp((17.269 * temp) / (237.3 + temp))
+    first_vel_corrn_ppm = param_c - (
+                (param_d * pressure) / (273.15 + temp) + ((11.27 * part_h2o_press) * (0.061078 * rel_humidity)) / (
+                    273.15 + temp))
+    first_vel_corrn_metres = dist * first_vel_corrn_ppm * (10 ** -6)
+    return first_vel_corrn_metres
+
+
 def va_conv(verta_hp, slope_dist, height_inst=0, height_tgt=0):
     """
     Function to convert vertical angles (zenith distances) and slope distances

geodepy/surveyconvert/classtools.py

@@ -8,6 +8,7 @@
 import itertools
 import operator
 from geodepy.convert import DMSAngle
+from geodepy.survey import first_vel_corrn
 
 
 class Coordinate(object):
@@ -232,3 +233,18 @@ def reducesetup(obslist, strict=False, zerodist=False):
     # Order list of meaned obs
     sorted_meanedobs = sorted(meanedobs, key=operator.attrgetter('hz_obs'))
     return sorted_meanedobs
+
+
+def first_vel_observations(obslist, params, temp, pressure, rel_humidity):
+    """
+    Performs a first velocity correction for all observed slope distances in a list of Observations
+    :param obslist: List of Observations (i.e. from one InstSetup)
+    :param params: Tuple of First Velocity Parameters C and D (see function first_vel_params)
+    :param temp: Observed Temperature (degrees Celsius)
+    :param pressure: Observed Pressure (hectopascals or millibars)
+    :param rel_humidity: Observed Relative Humidity (percentage)
+    :return: List of Observations with First Velocity Correction applied
+    """
+    for obs in obslist:
+        obs.sd_obs = obs.sd_obs + first_vel_corrn(obs.sd_obs, params, temp, pressure, rel_humidity)
+    return obslist

geodepy/surveyconvert/config.py

@@ -65,3 +65,17 @@ def removeobs(cfg_list, project):
                 elif remove_id == obs.to_id:
                     del setup.observation[num]
     return project
+
+
+def first_vel_cfg(cfg_list):
+    # Find First Velocity Parameter Inputs in cfg_list
+    for group in cfg_list:
+        group_header = group[0].lower()
+        if group_header.startswith('first vel'):
+            wavelength = float(group[1])
+            temperature = float(group[2])
+            pressure = float(group[3])
+            rel_humidity = float(group[4])
+            return wavelength, temperature, pressure, rel_humidity
+    else:
+        return None

geodepy/surveyconvert/fbk.py

@@ -8,8 +8,9 @@
 import os
 import numpy as np
 from geodepy.convert import DMSAngle
-from geodepy.surveyconvert.config import readconfig, renameobs, removeobs
-from geodepy.surveyconvert.classtools import Coordinate, InstSetup, Observation, reducesetup
+from geodepy.survey import first_vel_params
+from geodepy.surveyconvert.config import readconfig, renameobs, removeobs, first_vel_cfg
+from geodepy.surveyconvert.classtools import Coordinate, InstSetup, Observation, reducesetup, first_vel_observations
 from geodepy.surveyconvert.dna import dnaout_dirset, dnaout_va, dnaout_sd
 
 
@@ -32,10 +33,23 @@ def fbk2msr(path, cfg_path, strict=False, zerodist=False, same_stdev=False):
     fbk_project = renameobs(cfg, fbk_project)
     # Remove obs as per config file
     fbk_project = removeobs(cfg, fbk_project)
+    # Get First Velocity Correction Observations
+    first_vel_obs = first_vel_cfg(cfg)
     # Reduce observations in setups
     for setup in fbk_project:
         reduced_obs = reducesetup(setup.observation, strict, zerodist)
         setup.observation = reduced_obs
+    # Perform First Velocity Correction
+    if first_vel_obs is not None:
+        # Use wavelength and standard atmospheric parameters to get Parameters C and D
+        params = first_vel_params(first_vel_obs[0])
+        for setup in fbk_project:
+            corrected_obs = first_vel_observations(setup.observation,
+                                                   params,
+                                                   first_vel_obs[1],  # Observed Temperature
+                                                   first_vel_obs[2],  # Observed Pressure
+                                                   first_vel_obs[3])  # Observed Relative Humidity
+            setup.observation = corrected_obs
     # Produce Measurement format data from setups
     msr_raw = []
     for setup in fbk_project:

geodepy/surveyconvert/gsi.py

@@ -9,8 +9,9 @@
 import os
 from datetime import datetime
 from geodepy.convert import DMSAngle
-from geodepy.surveyconvert.config import readconfig, renameobs, removeobs
-from geodepy.surveyconvert.classtools import Coordinate, InstSetup, Observation, reducesetup
+from geodepy.survey import first_vel_params
+from geodepy.surveyconvert.config import readconfig, renameobs, removeobs, first_vel_cfg
+from geodepy.surveyconvert.classtools import Coordinate, InstSetup, Observation, reducesetup, first_vel_observations
 from geodepy.surveyconvert.dna import dnaout_dirset, dnaout_va, dnaout_sd
 
 
@@ -29,10 +30,25 @@ def gsi2msr(path, cfg_path=None):
         gsi_project = renameobs(cfg, gsi_project)
         # Remove obs as per config file
         gsi_project = removeobs(cfg, gsi_project)
+        # Get First Velocity Correction Observations
+        first_vel_obs = first_vel_cfg(cfg)
+    else:
+        first_vel_obs = None
     # Reduce observations in setups
     for setup in gsi_project:
         reduced_obs = reducesetup(setup.observation, strict=False, zerodist=False)
         setup.observation = reduced_obs
+    # Perform First Velocity Correction
+    if first_vel_obs is not None:
+        # Use wavelength and standard atmospheric parameters to get Parameters C and D
+        params = first_vel_params(first_vel_obs[0])
+        for setup in gsi_project:
+            corrected_obs = first_vel_observations(setup.observation,
+                                                   params,
+                                                   first_vel_obs[1],  # Observed Temperature
+                                                   first_vel_obs[2],  # Observed Pressure
+                                                   first_vel_obs[3])  # Observed Relative Humidity
+            setup.observation = corrected_obs
     # Produce Measurement format data from setups
     msr_raw = []
     for setup in gsi_project:
@@ -53,7 +69,7 @@ def gsi2msr(path, cfg_path=None):
             elif line.startswith('S'):
                 sdcount += 1
     obscount = dircount + vacount + sdcount
-    now = datetime.now()
+    now = datetime(2020, 1, 1)
     date = (str(now.day).rjust(2, '0') + '.'
             + str(now.month).rjust(2, '0') + '.'
             + str(now.year))

geodepy/tests/resources/ST0618HZ.msr

@@ -1,4 +1,4 @@
-!#=DNA 3.01 MSR    03.04.2019         GDA94    03.04.2019        173
+!#=DNA 3.01 MSR    01.01.2020         GDA94    01.01.2020        173
 D AU45                AU52                14                                 18 42 37.800  2.9335   
 D                                         AU52                               18 42 38.000  2.9335   
 D                                         AU52                               18 42 38.100  2.9335   

geodepy/tests/test_survey.py

@@ -0,0 +1,29 @@
+import unittest
+from geodepy.survey import first_vel_params, first_vel_corrn
+
+
+class TestSurveyConvert(unittest.TestCase):
+    def test_first_vel_params(self):
+        temperature = 12
+        pressure = 1013.25
+        relative_humidity = 60
+        edm_wavelength = 0.850
+        param_c = 281.781
+        param_d = 79.393
+        params_new = first_vel_params(edm_wavelength, temperature, pressure, relative_humidity)
+        self.assertEqual(round(params_new[0], 3), param_c)
+        self.assertEqual(round(params_new[1], 3), param_d)
+
+    def test_first_vel_corrn(self):
+        params = first_vel_params(0.85)
+        raw_obs_distance = 1117.8517
+        obs_temperature = 6.8
+        obs_pressure = 960.8
+        obs_relative_humidity = 58.6
+        correction = first_vel_corrn(raw_obs_distance, params, obs_temperature, obs_pressure, obs_relative_humidity)
+        corrected_obs_distance = raw_obs_distance + correction
+        self.assertEqual(round(corrected_obs_distance, 4), 1117.8618)
+
+
+if __name__ == '__main__':
+    unittest.main()

geodepy/tests/test_surveyconvert.py

@@ -53,4 +53,4 @@ def test_gsi2msr(self):
 
 
 if __name__ == '__main__':
-    unittest.main()
+    unittest.main()