Surveyconvert gsi2stn

geodepy/convert.py

@@ -5,7 +5,7 @@
 Conversions Module
 """
 
-from math import modf
+from math import modf, sin, cos, atan2, radians, degrees, sqrt
 
 
 def dec2hp(dec):
@@ -34,20 +34,26 @@ def __repr__(self):
     def __add__(self, other):
         return dec2dms(self.dec() + other.dec())
 
+    def __radd__(self, other):
+        return dec2dms(other.dec() + self.dec())
+
     def __sub__(self, other):
         return dec2dms(self.dec() - other.dec())
 
+    def __rsub__(self, other):
+        return dec2dms(other.dec() - self.dec())
+
     def __mul__(self, other):
         try:
             return dec2dms(self.dec() * other)
-        except ValueError:
-            raise ValueError('Multiply only defined between DMSAngle Object and Int or Float')
+        except TypeError:
+            raise TypeError('Multiply only defined between DMSAngle Object and Int or Float')
 
     def __rmul__(self, other):
         try:
             return dec2dms(other * self.dec())
-        except ValueError:
-            raise ValueError('Multiply only defined between DMSAngle Object and Int or Float')
+        except TypeError:
+            raise TypeError('Multiply only defined between DMSAngle Object and Int or Float')
 
     def __truediv__(self, other):
         return dec2dms(self.dec() / other)
@@ -103,14 +109,14 @@ def __sub__(self, other):
     def __mul__(self, other):
         try:
             return dec2ddm(self.dec() * other)
-        except ValueError:
-            raise ValueError('Multiply only defined between DMSAngle Object and Int or Float')
+        except TypeError:
+            raise TypeError('Multiply only defined between DMSAngle Object and Int or Float')
 
     def __rmul__(self, other):
         try:
             return dec2ddm(other * self.dec())
-        except ValueError:
-            raise ValueError('Multiply only defined between DMSAngle Object and Int or Float')
+        except TypeError:
+            raise TypeError('Multiply only defined between DMSAngle Object and Int or Float')
 
     def __truediv__(self, other):
         return dec2ddm(self.dec() / other)
@@ -127,6 +133,12 @@ def __eq__(self, other):
     def __ne__(self, other):
         return self.dec() != other.dec()
 
+    def __lt__(self, other):
+        return self.dec() < other.dec()
+
+    def __gt__(self, other):
+        return self.dec() > other.dec()
+
     def dec(self):
         if self.degree >= 0:
             return self.degree + (self.minute / 60)
@@ -171,6 +183,22 @@ def hp2ddm(hp):
     return DDMAngle(degree, minute) if hp >= 0 else DDMAngle(-degree, minute)
 
 
+def polar2rect(r, theta):
+    x = r * sin(radians(theta))
+    y = r * cos(radians(theta))
+    return x, y
+
+
+def rect2polar(x, y):
+    r = sqrt(x ** 2 + y ** 2)
+    theta = atan2(x, y)
+    if theta < 0:
+        theta = degrees(theta) + 360
+    else:
+        theta = degrees(theta)
+    return r, theta
+
+
 # ----------------
 # Legacy Functions
 # dd2dms refactored as dec2hp

geodepy/survey.py

@@ -5,8 +5,9 @@
 Survey Module
 """
 
-from math import sqrt, sin, cos, asin, radians, degrees, exp
-from geodepy.convert import hp2dec, dec2hp
+from math import sqrt, sin, cos, atan, radians, degrees, exp
+from statistics import mean, stdev
+from geodepy.convert import rect2polar, polar2rect
 
 
 def first_vel_params(wavelength, temp=12, pressure=1013.25, rel_humidity=60):
@@ -49,48 +50,77 @@ def first_vel_corrn(dist, first_vel_param, temp, pressure, rel_humidity):
     return first_vel_corrn_metres
 
 
-def va_conv(verta_hp, slope_dist, height_inst=0, height_tgt=0):
+def precise_inst_ht(vert_list, spacing, offset):
+    """
+    Uses a set of Vertical Angle Observations taken to a
+    levelling staff at regular intervals to determine the
+    height of the instrument above a reference mark
+    :param vert_list: List of Vertical (Zenith) Angle Observations (minimum of 3) in Decimal Degrees format
+    :param spacing: Distance in metres between each vertical angle observation
+    :param offset: Lowest observed height above reference mark
+    :return: Instrument Height above reference mark and its standard deviation
+    """
+    if len(vert_list) < 3:
+        raise ValueError('ValueError: 3 or more vertical angles required')
+    vert_list.sort(reverse=True)
+    vert_pairs = [(va1, va2) for va1, va2 in zip(vert_list, vert_list[1:])]
+    base_ht = []
+    height_comp = []
+    for num, pair in enumerate(vert_pairs):
+        base_ht_pair = offset + num * spacing
+        base_ht.append(base_ht_pair)
+        dist_a = sin(radians(pair[1])) * (spacing / (sin(radians(pair[0] - pair[1]))))
+        delta_ht = dist_a * (sin(radians(pair[0] - 90)))
+        height_comp.append(delta_ht + base_ht[num])
+    return round(mean(height_comp), 5), round(stdev(height_comp), 5)
+
+
+def joins(east1, north1, east2, north2):
+    return rect2polar(east2 - east1, north2 - north1)
+
+
+def radiations(east1, north1, brg1to2, dist, rotation=0, psf=1):
+    delta_east, delta_north = polar2rect(dist * psf, brg1to2 + rotation)
+    return east1 + delta_east, north1 + delta_north
+
+
+def va_conv(zenith_angle, slope_dist, height_inst=0, height_tgt=0):
     """
     Function to convert vertical angles (zenith distances) and slope distances
     into horizontal distances and changes in height. Instrument and Target
     heights can be entered to allow computation of zenith and slope distances
     between ground points.
 
-    :param verta_hp:        Vertical Angle from Instrument to Target, expressed
-                            in HP Format (DDD.MMSSSSSS)
+    :param zenith_angle:    Zenith Angle from Instrument to Target, expressed
+                            in decimal degrees
     :param slope_dist:      Slope Distance from Instrument to Target in metres
     :param height_inst:     Height of Instrument. Optional - Default Value of 0m
     :param height_tgt:      Height of Target. Optional - Default Value of 0m
 
-    :return: verta_pt_hp:   Vertical Angle between Ground Points, expressed
-                            in HP Format (DDD.MMSSSSSS)
+    :return: vert_angle_pt: Vertical Angle between Ground Points, expressed
+                            in decimal degrees
     :return: slope_dist_pt: Slope Distance between Ground Points in metres
     :return: hz_dist:       Horizontal Distance
     :return: delta_ht:      Change in height between Ground Points in metres
     """
     # Convert Zenith Angle to Vertical Angle
     try:
-        if verta_hp == 0 or verta_hp == 180:
+        if zenith_angle == 0 or zenith_angle == 180:
             raise ValueError
-        elif 0 < verta_hp < 180:
-            verta = radians(90 - hp2dec(verta_hp))
-        elif 180 < verta_hp < 360:
-            verta = radians(270 - hp2dec(verta_hp))
+        elif 0 < zenith_angle < 180:
+            zenith_angle = radians(90 - zenith_angle)
+        elif 180 < zenith_angle < 360:
+            zenith_angle = radians(270 - zenith_angle)
         else:
             raise ValueError
     except ValueError:
-        print('ValueError: Vertical Angle Invalid')
-        return
+        raise ValueError('ValueError: Vertical Angle Invalid')
     # Calculate Horizontal Dist and Delta Height
-    hz_dist = slope_dist * cos(verta)
-    delta_ht = slope_dist * sin(verta)
+    hz_dist = slope_dist * cos(zenith_angle)
+    delta_ht = slope_dist * sin(zenith_angle)
     # Account for Target and Instrument Heights
-    if height_inst == 0 and height_tgt == 0:
-        verta_pt_hp = verta_hp
-        slope_dist_pt = slope_dist
-    else:
-        delta_ht = height_inst + delta_ht - height_tgt
-        slope_dist_pt = sqrt(delta_ht ** 2 + hz_dist ** 2)
-        verta_pt = asin(delta_ht / slope_dist)
-        verta_pt_hp = dec2hp(degrees(verta_pt) + 90)
-    return verta_pt_hp, slope_dist_pt, hz_dist, delta_ht
+    delta_ht = height_inst + delta_ht - height_tgt
+    slope_dist_pt = sqrt(delta_ht ** 2 + hz_dist ** 2)
+    vert_angle_pt = atan(delta_ht / hz_dist)
+    vert_angle_pt = degrees(vert_angle_pt)
+    return vert_angle_pt, slope_dist_pt, hz_dist, delta_ht

geodepy/surveyconvert/classtools.py

@@ -7,7 +7,8 @@
 
 import itertools
 import operator
-from geodepy.convert import DMSAngle
+from statistics import mean
+from geodepy.convert import DMSAngle, dec2dms
 from geodepy.survey import first_vel_corrn
 
 
@@ -55,13 +56,14 @@ def __iter__(self):
 
 class Observation(object):
     def __init__(self, from_id, to_id, inst_height=0.0, target_height=0.0,
-                 face='FL', hz_obs=DMSAngle(0, 0, 0), va_obs=DMSAngle(0, 0, 0),
+                 face='FL', rounds=0.5, hz_obs=DMSAngle(0, 0, 0), va_obs=DMSAngle(0, 0, 0),
                  sd_obs=0.0, hz_dist=0.0, vert_dist=0.0):
         self.from_id = from_id
         self.to_id = to_id
         self.inst_height = inst_height
         self.target_height = target_height
         self.face = face
+        self.rounds = rounds
         self.hz_obs = hz_obs
         self.va_obs = va_obs
         self.sd_obs = sd_obs
@@ -74,6 +76,7 @@ def __repr__(self):
                 + '; inst_height ' + repr(self.inst_height)
                 + '; target_height ' + repr(self.target_height)
                 + '; face ' + repr(self.face)
+                + '; rounds ' + repr(self.rounds)
                 + '; hz_obs ' + repr(self.hz_obs)
                 + '; va_obs ' + repr(self.va_obs)
                 + '; sd_obs ' + repr(self.sd_obs)
@@ -103,12 +106,14 @@ def changeface(self):
                            self.inst_height,
                            self.target_height,
                            newface,
+                           self.rounds,
                            hz_switch,
                            va_switch,
                            self.sd_obs,
                            self.hz_obs,
                            self.vert_dist)
 
+
 def meanfaces(ob1, ob2):
     """
     Take two Observations and return their mean Face Left Sense Observation
@@ -172,12 +177,13 @@ def meanfaces(ob1, ob2):
                            ob1.inst_height,
                            ob1.target_height,
                            ob1.face,
+                           ob1.rounds + ob2.rounds,
                            meaned_hz,
                            meaned_va,
                            meaned_sd)
 
 
-def reducesetup(obslist, strict=False, zerodist=False):
+def reducesetup(obslist, strict=False, zerodist=False, meanmulti=False):
     """
     Takes a list of Observations from one setup and
     means together FL, FR pairs of Observations.
@@ -186,9 +192,10 @@ def reducesetup(obslist, strict=False, zerodist=False):
     single-face Obs are included and converted to Face Left
     :param zerodist: If True, Obs with Slope Distance of Zero are included.
     If False, these are ignored
+    :param meanmulti: If True, multiple rounds of observations are reduced
+    to a single FL-sense Observation
     :return: a reduced list of Observations
     """
-
     # Remove obs with sd_obs == 0
     if not zerodist:
         for ob in obslist:
@@ -230,6 +237,50 @@ def reducesetup(obslist, strict=False, zerodist=False):
                 for pair in pairedlist:
                     meanob = meanfaces(pair[0], pair[1])
                     meanedobs.append(meanob)
+    # Mean multiple repeated measurements into a single FL observation
+    if meanmulti:
+        multiob = []
+        to_ids = set([i.to_id for i in meanedobs])
+        for id in to_ids:
+            matchedobs = []
+            for ob in meanedobs:
+                if ob.to_id == id:
+                    matchedobs.append(ob)
+            ob1 = matchedobs[0]
+            repeat_hz = [ob.hz_obs.dec() for ob in matchedobs]
+            repeat_hz = sorted(repeat_hz)
+            repeat_va = [ob.va_obs.dec() for ob in matchedobs]
+            repeat_sd = [ob.sd_obs for ob in matchedobs]
+            # Finds where Horizontal Obseravtions are either side of 360, converts those on 360 side to negative
+            hz_diff = [j - i for i, j in zip(repeat_hz[:-1], repeat_hz[1:])]
+            bigjump = 1e100
+            for num, i in enumerate(hz_diff):
+                if i > 180:
+                    bigjump = num
+            for num, ob in enumerate(repeat_hz):
+                if num > bigjump:
+                    repeat_hz[num] = repeat_hz[num] - 360
+            # Mean Repeated Observations
+            mean_hz = mean(repeat_hz)
+            if mean_hz < 0:
+                mean_hz + 360
+            mean_va = mean(repeat_va)
+            mean_sd = mean(repeat_sd)
+            # Compute number of rounds of observations completed
+            sum_rounds = 0
+            for ob in matchedobs:
+                sum_rounds += ob.rounds
+            # Output Meaned Observation
+            multiob.append(Observation(ob1.from_id,
+                                       id,
+                                       ob1.inst_height,
+                                       ob1.target_height,
+                                       ob1.face,
+                                       sum_rounds,
+                                       dec2dms(mean_hz),
+                                       dec2dms(mean_va),
+                                       mean_sd))
+        meanedobs = multiob
     # Order list of meaned obs
     sorted_meanedobs = sorted(meanedobs, key=operator.attrgetter('hz_obs'))
     return sorted_meanedobs

geodepy/surveyconvert/config.py

@@ -79,3 +79,17 @@ def first_vel_cfg(cfg_list):
             return wavelength, temperature, pressure, rel_humidity
     else:
         return None
+
+
+def stdev_cfg(cfg_list):
+    # Find Survey Apriori Standard Deviations in cfg_list
+    for group in cfg_list:
+        group_header = group[0].lower()
+        if group_header.startswith('standard dev'):
+            stdev_angular = float(group[1])
+            stdev_pointing = float(group[2])
+            stdev_distconst = float(group[3])
+            stdev_distppm = float(group[4])
+            return stdev_angular, stdev_pointing, stdev_distconst, stdev_distppm
+    else:
+        return None