From 7b4aa150a3731ecf51e3963466e1ed35994c744d Mon Sep 17 00:00:00 2001
From: pcain-usgs <pcain@usgs.gov>
Date: Thu, 9 Apr 2020 12:27:06 -0600
Subject: [PATCH] Clean up/ add comments
---
geomagio/residual/Calculation.py | 176 +++++++++---------
geomagio/residual/Reading.py | 2 +-
.../residual/SpreadsheetAbsolutesFactory.py | 6 +-
3 files changed, 91 insertions(+), 93 deletions(-)
diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py
index 4d407f5a6..86f3cc954 100644
--- a/geomagio/residual/Calculation.py
+++ b/geomagio/residual/Calculation.py
@@ -22,10 +22,7 @@ class Calculate(object):
angle: float = None,
residual: float = None,
ordinate: Ordinate = None,
- baseline: float = None,
f: float = None,
- inclination: float = None,
- hs: int = None,
ud: int = None,
shift: int = None,
pm: int = None,
@@ -33,26 +30,26 @@ class Calculate(object):
self.angle = angle
self.residual = residual
self.ordinate = ordinate
- self.baseline = baseline
self.f = f
- self.inclination = inclination
- self.hs = hs
self.ud = ud
self.shift = shift
self.pm = pm
def calculate(Reading):
- # get average ordinate values across h, e, z, and f
+ """
+ Calculate/recalculate absolute from a Reading object's
+ ordinates, measurements, and metadata.
+ Outputs a list of absolutes containing baseline, absolute,
+ and element name. Also reutrns the scale value.
+ """
+ # define measurement types used to calculate inclination
+ inclination_types = [mt.NORTH_DOWN, mt.NORTH_UP, mt.SOUTH_DOWN, mt.SOUTH_UP]
+ # get ordinate values across h, e, z, and f for inclination measurement types
inclination_ordinates = [
- o
- for o in Reading.ordinates
- if "West" not in o.measurement_type.capitalize()
- and "East" not in o.measurement_type.capitalize()
- and "NorthDownScale" not in o.measurement_type.capitalize()
+ o for o in Reading.ordinates if o.measurement_type in inclination_types
]
-
- inclination_ordinates = inclination_ordinates[0:-2]
+ # get average ordinate values across h, e, z, and f
mean = average_ordinate(inclination_ordinates, None)
# calculate inclination
inclination, f = calculate_I(
@@ -63,13 +60,10 @@ def calculate(Reading):
Reading.metadata,
)
# calculate absolutes
- # FIXME: change to self.pier_correction
- Habs, Zabs = calculate_absolutes(
- f, inclination, Reading.metadata["pier_correction"]
- )
- # calculate baselines
+ Habs, Zabs = calculate_absolutes(f, inclination)
+ # calculate baselines for H and Z
Hb, Zb = calculate_baselines(Habs, Zabs, mean, Reading.pier_correction)
- # calculate scale value for declination
+ # calculate scale value
scale_ordinates = Reading.ordinate_index()["NorthDownScale"]
scale_measurements = Reading.measurement_index()["NorthDownScale"]
scale = calculate_scale(
@@ -77,9 +71,8 @@ def calculate(Reading):
scale_ordinates,
scale_measurements,
inclination,
- Reading.metadata["pier_correction"],
)
- # calculate declination and
+ # calculate declination absolute and baseline
Db, Dabs = calculate_D(
Reading.ordinate_index(),
Reading.measurements,
@@ -104,11 +97,10 @@ def calculate_I(measurements, ordinates, ordinates_index, mean, metadata):
average ordinates from every measurement, and metadata.
Returns inclination angle and calculated average f
"""
- # get first inclination angle, assumed to be southdown
+ # get first inclination angle, assumed to be the southdown angle
Iprime = average_angle(measurements, "SouthDown")
if Iprime >= 100:
Iprime -= 200
- # Iprime = (np.pi / 200)*(Iprime)
# get multiplier for hempisphere the observatory is located in
# 1 if observatory is in northern hemisphere
# -1 if observatory is in southern hemisphere
@@ -117,7 +109,6 @@ def calculate_I(measurements, ordinates, ordinates_index, mean, metadata):
southdown = Calculate(
shift=-200,
ud=1,
- hs=hs,
pm=1,
angle=average_angle(measurements, "SouthDown"),
residual=average_residual(measurements, "SouthDown"),
@@ -127,7 +118,6 @@ def calculate_I(measurements, ordinates, ordinates_index, mean, metadata):
southup = Calculate(
shift=200,
ud=-1,
- hs=hs,
pm=-1,
angle=average_angle(measurements, "SouthUp"),
residual=average_residual(measurements, "SouthUp"),
@@ -137,7 +127,6 @@ def calculate_I(measurements, ordinates, ordinates_index, mean, metadata):
northup = Calculate(
shift=0,
ud=-1,
- hs=hs,
pm=1,
angle=average_angle(measurements, "NorthUp"),
residual=average_residual(measurements, "NorthUp"),
@@ -147,33 +136,33 @@ def calculate_I(measurements, ordinates, ordinates_index, mean, metadata):
northdown = Calculate(
shift=400,
ud=1,
- hs=hs,
pm=-1,
angle=average_angle(measurements, "NorthDown"),
residual=average_residual(measurements, "NorthDown"),
ordinate=average_ordinate(ordinates_index, "NorthDown"),
)
-
+ # set inclination value for looping = Iprime
inclination = Iprime
+ # add one to inclination value to enter the loop
Inclination = inclination + 1
-
+ # loop condition
while abs(Inclination - inclination) > 0.0001:
+ # set temporary inclination variable to hold previous step's inclination
Inclination = inclination
-
+ # calculate f for inclination measurement types
southdown.f = calculate_f(southdown.ordinate, mean, inclination)
northdown.f = calculate_f(northdown.ordinate, mean, inclination)
southup.f = calculate_f(southup.ordinate, mean, inclination)
northup.f = calculate_f(northup.ordinate, mean, inclination)
-
- northup.inclination = calculate_measurement_inclination(northup)
- southdown.inclination = calculate_measurement_inclination(southdown)
- northdown.inclination = calculate_measurement_inclination(northdown)
- southup.inclination = calculate_measurement_inclination(southup)
-
+ # gather measurements into array(necessary because f were recalculated)
measurements = [southup, southdown, northup, northdown]
-
- inclination = np.average([i.inclination for i in measurements])
+ # average f values for inclination measurement types
f = np.average([i.f for i in measurements])
+ # calculation inclination for each inclination measurement type and average
+ inclination = np.average(
+ [calculate_measurement_inclination(m, hs) for m in measurements]
+ )
+ # loop exits once the difference in average inclination between steps is lower than 0.0001
return inclination, f
@@ -184,15 +173,24 @@ def calculate_D(ordinates_index, measurements, measurements_index, AZ, Hb):
ordinates, measurements, measurement_index(dictionary), azimuth and H baseline
Returns absolute and baseline for declination.
"""
- # compute average angle from marks
+ # specify mark measurement types
+ mark_types = [
+ mt.FIRST_MARK_DOWN,
+ mt.FIRST_MARK_UP,
+ mt.SECOND_MARK_DOWN,
+ mt.SECOND_MARK_UP,
+ ]
+ # average mark angles
+ # note that angles are being converted to geon
average_mark = np.average(
[
convert_to_geon(m.angle)
for m in measurements
- if "mark" in m.measurement_type.capitalize()
+ if m.measurement_type in mark_types
]
)
-
+ # add 100 if mark up is greater than mark down
+ # subtract 100 otherwise
if (
measurements_index["FirstMarkUp"][0].angle
< measurements_index["FirstMarkDown"][0].angle
@@ -201,67 +199,66 @@ def calculate_D(ordinates_index, measurements, measurements_index, AZ, Hb):
else:
average_mark -= 100
- # if average_mark < 200:
- # average_mark += 100
- # else:
- # average_mark -= 100
-
- # gather calculation objects for each measurement type
+ # gather calculation objects for each declination measurement type
+ # note that the pm(plus minus) multiplier has been repurposed.
+ # West facing measurements have a multiplier of -1
+ # East facing measurements have a multipllier of 1
westdown = Calculate(
- baseline=Hb,
angle=average_angle(measurements_index, "WestDown"),
residual=average_residual(measurements_index, "WestDown"),
ordinate=average_ordinate(ordinates_index, "WestDown"),
- ud=-1,
+ pm=-1,
)
westup = Calculate(
- baseline=Hb,
angle=average_angle(measurements_index, "WestUp"),
residual=average_residual(measurements_index, "WestUp"),
ordinate=average_ordinate(ordinates_index, "WestUp"),
- ud=-1,
+ pm=-1,
)
eastdown = Calculate(
- baseline=Hb,
angle=average_angle(measurements_index, "EastDown"),
residual=average_residual(measurements_index, "EastDown"),
ordinate=average_ordinate(ordinates_index, "EastDown"),
- ud=1,
+ pm=1,
)
eastup = Calculate(
- baseline=Hb,
angle=average_angle(measurements_index, "EastUp"),
residual=average_residual(measurements_index, "EastUp"),
ordinate=average_ordinate(ordinates_index, "EastUp"),
- ud=1,
+ pm=1,
)
- # AZ = convert_to_geon(AZ, include_seconds=False)
+ # convert azimuth to geon
AZ = (int(AZ / 100) + (AZ % 100) / 60) / 0.9
- # gather measurements into array
+ # gather declination measurements into array
measurements = [westdown, westup, eastdown, eastup]
- # get average meridian angle from measurement types
- meridian = np.average([calculate_meridian_term(i) for i in measurements])
- # add average mark, meridian, and azimuth angle to get declination baseline
- Db_abs = (meridian - average_mark) + AZ
- Db = round(Db_abs * 54, 2)
+ # average meridian terms calculated from each declination measurements
+ meridian = np.average([calculate_meridian_term(i, Hb) for i in measurements])
+ # add subtract average mark angle from average meridian angle and add
+ # azimuth(in geon) to get the declination baseline
+ D_b = (meridian - average_mark) + AZ
+ # convert decimal baseline into dms baseline
+ Db = round(D_b * 54, 2)
Db_dms = int(Db / 60) * 100 + ((Db / 60) % 1) * 60
- # calculate declination absolute
+ # gather first declination measurements' H ans E ordinates
wd_E_1 = ordinates_index["WestDown"][0].e
wd_H_1 = ordinates_index["WestDown"][0].h
- Dabs = Db_abs + np.arctan(wd_E_1 / (Hb + wd_H_1)) * (200 / np.pi)
+ # calculate Declination baseline
+ Dabs = D_b + np.arctan(wd_E_1 / (Hb + wd_H_1)) * (200 / np.pi)
Dabs = round(Dabs * 54, 1)
+ # convert decimal absolute into dms absolute
Dabs_dms = int(Dabs / 60) * 100 + ((Dabs / 60) % 1) * 60
return Db_dms, Dabs_dms
-def calculate_absolutes(f, inclination, pier_correction):
+def calculate_absolutes(f, inclination):
"""
Calculate absolutes for H, Z and F from computed
- average f value(from inclination computations),
- calculated inclination angle, and pier correction(metadata).
- Returns baselines for H, Z, and F
+ average f value(from inclination computations) and
+ calculated inclination angle.
+ Returns baselines for H and Z
"""
+ # convert inclination to radians
i = (np.pi / 200) * (inclination)
Habs = f * np.cos(i)
Zabs = f * np.sin(i)
@@ -271,11 +268,11 @@ def calculate_absolutes(f, inclination, pier_correction):
def calculate_baselines(Habs, Zabs, mean, pier_correction):
"""
- Calculate baselines with H and Z absolutes, and
- average ordinates across all measurements.
+ Calculate baselines with H and Z absolutes,
+ average ordinates across all measurements,
+ and pier correction(metadata).
Returns H and Z baselines
"""
- # FIXME: Figure out where 0.3 comes from
correction = pier_correction / 5
Hb = round(np.sqrt(Habs ** 2 - mean.e ** 2) - mean.h, 1) - correction
Zb = round(Zabs - mean.z, 1) - correction
@@ -283,10 +280,12 @@ def calculate_baselines(Habs, Zabs, mean, pier_correction):
return Hb, Zb
-def calculate_scale(f, ordinates, measurements, inclination, pier_correction):
+def calculate_scale(f, ordinates, measurements, inclination):
"""
Calculate scale value from calulated f(from inclination computations),
- calculated inclination, and pier correction(metadata)
+ inclination, and the measurements/ordinates taken for scaling purposes.
+ Such measurements usually present themselves as a set of three North Down
+ measurements, where the final two measuremets(and ordinates) are used for scaling.
"""
first_ord = ordinates[0]
second_ord = ordinates[1]
@@ -316,7 +315,6 @@ def average_angle(measurements, type):
Compute average angle from a dictionary of
measurements and specified measurement type.
"""
- # FIXME: change repetitive checks
return np.average([convert_to_geon(m.angle) for m in measurements[type]])
@@ -333,7 +331,6 @@ def average_ordinate(ordinates, type):
Compute average ordinate from a dictionary
of ordinates and specified measurement type.
"""
- # FIXME: change repetitive checks
if type is not None:
ordinates = ordinates[type]
o = Ordinate(measurement_type=type)
@@ -349,7 +346,6 @@ def calculate_f(ordinate, mean, inclination):
and calculated inclination.
"""
# get channel means form all ordinates
- # FIXME: don't unpack ordinates
# calculate f using current step's inclination angle
f = (
mean.f
@@ -360,39 +356,38 @@ def calculate_f(ordinate, mean, inclination):
return f
-def calculate_measurement_inclination(calc):
+def calculate_measurement_inclination(calculation, hs):
"""
Calculate a measurement's inclination value using
Calculate items' elements.
"""
- return calc.shift + calc.pm * (
- +calc.angle
- + calc.ud * (calc.hs * np.arcsin(calc.residual / calc.f) * 200 / np.pi)
+ return calculation.shift + calculation.pm * (
+ +calculation.angle
+ + calculation.ud
+ * (hs * np.arcsin(calculation.residual / calculation.f) * 200 / np.pi)
)
-def calculate_meridian_term(calculation):
+def calculate_meridian_term(calculation, Hb):
"""
Calculate meridian value from a measurement type
using a Calculate object and H's baseline value.
"""
A1 = np.arcsin(
calculation.residual
- / np.sqrt(
- (calculation.ordinate.h + calculation.baseline) ** 2
- + (calculation.ordinate.e) ** 2
- )
- )
- A2 = np.arctan(
- calculation.ordinate.e / (calculation.ordinate.h + calculation.baseline)
+ / np.sqrt((calculation.ordinate.h + Hb) ** 2 + (calculation.ordinate.e) ** 2)
)
+ A2 = np.arctan(calculation.ordinate.e / (calculation.ordinate.h + Hb))
A1 = (200 / np.pi) * (A1)
A2 = (200 / np.pi) * (A2)
- meridian_term = calculation.angle + (calculation.ud * A1) - A2
+ meridian_term = calculation.angle + (calculation.pm * A1) - A2
return meridian_term
def convert_to_geon(angle, include_seconds=True):
+ """
+ Convert angles from measurements to geon
+ """
degrees = int(angle)
minutes = int((angle % 1) * 100) / 60
if include_seconds:
@@ -400,4 +395,5 @@ def convert_to_geon(angle, include_seconds=True):
else:
seconds = 0
dms = (degrees + minutes + seconds) / 0.9
+
return dms
diff --git a/geomagio/residual/Reading.py b/geomagio/residual/Reading.py
index aa4035114..531eb6a46 100644
--- a/geomagio/residual/Reading.py
+++ b/geomagio/residual/Reading.py
@@ -25,7 +25,7 @@ class Reading(BaseModel):
absolutes: Optional[List[Absolute]] = None
azimuth: float = 0
- hemisphere: float = 1 # maybe hemisphere should be calculated from latitude
+ hemisphere: float = 1
measurements: Optional[List[Measurement]] = []
ordinates: Optional[List[Ordinate]] = []
metadata: Optional[Dict] = []
diff --git a/geomagio/residual/SpreadsheetAbsolutesFactory.py b/geomagio/residual/SpreadsheetAbsolutesFactory.py
index 0fd164794..1952e8429 100644
--- a/geomagio/residual/SpreadsheetAbsolutesFactory.py
+++ b/geomagio/residual/SpreadsheetAbsolutesFactory.py
@@ -279,8 +279,10 @@ class SpreadsheetAbsolutesFactory(object):
and self._parse_measurements(measurement_sheet, metadata["date"],)
or None
)
- ordinates = include_measurements and self._parse_ordinates(
- measurement_sheet, metadata["date"]
+ ordinates = (
+ include_measurements
+ and self._parse_ordinates(measurement_sheet, metadata["date"])
+ or None
)
return Reading(
absolutes=absolutes,
--
GitLab