From 8f0cfdcb10528b72eb33bd70d9999f1cfcdcc70d Mon Sep 17 00:00:00 2001
From: Jeremy Fee <jmfee@usgs.gov>
Date: Tue, 28 Apr 2020 22:39:42 -0600
Subject: [PATCH] Additional refactor of calculations
---
Pipfile | 4 +-
Pipfile.lock | 163 ++++++++------
geomagio/residual/Calculation.py | 318 ++++++++++++++++-----------
geomagio/residual/Measurement.py | 21 +-
geomagio/residual/MeasurementType.py | 51 ++++-
geomagio/residual/Reading.py | 23 +-
geomagio/residual/__init__.py | 29 ++-
test/residual_test/residual_test.py | 7 +-
8 files changed, 378 insertions(+), 238 deletions(-)
diff --git a/Pipfile b/Pipfile
index 0b6c79757..36de01836 100644
--- a/Pipfile
+++ b/Pipfile
@@ -17,13 +17,15 @@ obspy = ">1.2.0"
pycurl = "*"
authlib = "*"
cryptography = "*"
-databases = {extras = ["sqlite"],version = "*"}
+databases = {extras = ["postgresql", "sqlite"], version = "*"}
fastapi = "*"
httpx = "==0.11.1"
openpyxl = "*"
pydantic = "==1.4"
sqlalchemy = "*"
sqlalchemy-utc = "*"
+typing = "*"
+typing-extensions = "*"
uvicorn = "*"
[pipenv]
diff --git a/Pipfile.lock b/Pipfile.lock
index c93459f43..a8cdfa11a 100644
--- a/Pipfile.lock
+++ b/Pipfile.lock
@@ -1,7 +1,7 @@
{
"_meta": {
"hash": {
- "sha256": "feb945aefc3117570ac55b98510d96243c861469466ee798fa8b49416693f89d"
+ "sha256": "809e5e04e4d1957e08e475b75c4c425c1b63b7d1fdaee30e79b077e84e115504"
},
"pipfile-spec": 6,
"requires": {},
@@ -103,35 +103,35 @@
},
"click": {
"hashes": [
- "sha256:8a18b4ea89d8820c5d0c7da8a64b2c324b4dabb695804dbfea19b9be9d88c0cc",
- "sha256:e345d143d80bf5ee7534056164e5e112ea5e22716bbb1ce727941f4c8b471b9a"
+ "sha256:d2b5255c7c6349bc1bd1e59e08cd12acbbd63ce649f2588755783aa94dfb6b1a",
+ "sha256:dacca89f4bfadd5de3d7489b7c8a566eee0d3676333fbb50030263894c38c0dc"
],
- "version": "==7.1.1"
+ "version": "==7.1.2"
},
"cryptography": {
"hashes": [
- "sha256:0cacd3ef5c604b8e5f59bf2582c076c98a37fe206b31430d0cd08138aff0986e",
- "sha256:192ca04a36852a994ef21df13cca4d822adbbdc9d5009c0f96f1d2929e375d4f",
- "sha256:19ae795137682a9778892fb4390c07811828b173741bce91e30f899424b3934d",
- "sha256:1b9b535d6b55936a79dbe4990b64bb16048f48747c76c29713fea8c50eca2acf",
- "sha256:2a2ad24d43398d89f92209289f15265107928f22a8d10385f70def7a698d6a02",
- "sha256:3be7a5722d5bfe69894d3f7bbed15547b17619f3a88a318aab2e37f457524164",
- "sha256:49870684da168b90110bbaf86140d4681032c5e6a2461adc7afdd93be5634216",
- "sha256:587f98ce27ac4547177a0c6fe0986b8736058daffe9160dcf5f1bd411b7fbaa1",
- "sha256:5aca6f00b2f42546b9bdf11a69f248d1881212ce5b9e2618b04935b87f6f82a1",
- "sha256:6b744039b55988519cc183149cceb573189b3e46e16ccf6f8c46798bb767c9dc",
- "sha256:6b91cab3841b4c7cb70e4db1697c69f036c8bc0a253edc0baa6783154f1301e4",
- "sha256:7598974f6879a338c785c513e7c5a4329fbc58b9f6b9a6305035fca5b1076552",
- "sha256:7a279f33a081d436e90e91d1a7c338553c04e464de1c9302311a5e7e4b746088",
- "sha256:95e1296e0157361fe2f5f0ed307fd31f94b0ca13372e3673fa95095a627636a1",
- "sha256:9fc9da390e98cb6975eadf251b6e5fa088820141061bf041cd5c72deba1dc526",
- "sha256:cc20316e3f5a6b582fc3b029d8dc03aabeb645acfcb7fc1d9848841a33265748",
- "sha256:d1bf5a1a0d60c7f9a78e448adcb99aa101f3f9588b16708044638881be15d6bc",
- "sha256:ed1d0760c7e46436ec90834d6f10477ff09475c692ed1695329d324b2c5cd547",
- "sha256:ef9a55013676907df6c9d7dd943eb1770d014f68beaa7e73250fb43c759f4585"
+ "sha256:091d31c42f444c6f519485ed528d8b451d1a0c7bf30e8ca583a0cac44b8a0df6",
+ "sha256:18452582a3c85b96014b45686af264563e3e5d99d226589f057ace56196ec78b",
+ "sha256:1dfa985f62b137909496e7fc182dac687206d8d089dd03eaeb28ae16eec8e7d5",
+ "sha256:1e4014639d3d73fbc5ceff206049c5a9a849cefd106a49fa7aaaa25cc0ce35cf",
+ "sha256:22e91636a51170df0ae4dcbd250d318fd28c9f491c4e50b625a49964b24fe46e",
+ "sha256:3b3eba865ea2754738616f87292b7f29448aec342a7c720956f8083d252bf28b",
+ "sha256:651448cd2e3a6bc2bb76c3663785133c40d5e1a8c1a9c5429e4354201c6024ae",
+ "sha256:726086c17f94747cedbee6efa77e99ae170caebeb1116353c6cf0ab67ea6829b",
+ "sha256:844a76bc04472e5135b909da6aed84360f522ff5dfa47f93e3dd2a0b84a89fa0",
+ "sha256:88c881dd5a147e08d1bdcf2315c04972381d026cdb803325c03fe2b4a8ed858b",
+ "sha256:96c080ae7118c10fcbe6229ab43eb8b090fccd31a09ef55f83f690d1ef619a1d",
+ "sha256:a0c30272fb4ddda5f5ffc1089d7405b7a71b0b0f51993cb4e5dbb4590b2fc229",
+ "sha256:bb1f0281887d89617b4c68e8db9a2c42b9efebf2702a3c5bf70599421a8623e3",
+ "sha256:c447cf087cf2dbddc1add6987bbe2f767ed5317adb2d08af940db517dd704365",
+ "sha256:c4fd17d92e9d55b84707f4fd09992081ba872d1a0c610c109c18e062e06a2e55",
+ "sha256:d0d5aeaedd29be304848f1c5059074a740fa9f6f26b84c5b63e8b29e73dfc270",
+ "sha256:daf54a4b07d67ad437ff239c8a4080cfd1cc7213df57d33c97de7b4738048d5e",
+ "sha256:e993468c859d084d5579e2ebee101de8f5a27ce8e2159959b6673b418fd8c785",
+ "sha256:f118a95c7480f5be0df8afeb9a11bd199aa20afab7a96bcf20409b411a3a85f0"
],
"index": "pypi",
- "version": "==2.9"
+ "version": "==2.9.2"
},
"cycler": {
"hashes": [
@@ -146,10 +146,10 @@
"sqlite"
],
"hashes": [
- "sha256:a04db1d158a91db7bd49db16e14266e8e6c7336f06f88c700147690683c769a3"
+ "sha256:92b7f66c74d21eba54899ba88352311019131fb01ceaf076452949429e8d1b12"
],
"index": "pypi",
- "version": "==0.2.6"
+ "version": "==0.3.1"
},
"decorator": {
"hashes": [
@@ -201,10 +201,10 @@
},
"hstspreload": {
"hashes": [
- "sha256:acbf4d6fd362b363ce567db56a8667c4f0e43073001add9c4e43f449c11e1d81",
- "sha256:ca0af58f803598a737d5f325212e488b24de1ac9f0cf6fa12da1a3f4651914e6"
+ "sha256:5e42c9c9b925d793d1732e20f59da7c9d6d340bde82ead7d49031bde59991219",
+ "sha256:f2571debd824aa5ec34816a2e08badda043181e908710317b9826b5c30d3bc17"
],
- "version": "==2020.4.14"
+ "version": "==2020.4.28"
},
"httptools": {
"hashes": [
@@ -327,30 +327,30 @@
},
"numpy": {
"hashes": [
- "sha256:1598a6de323508cfeed6b7cd6c4efb43324f4692e20d1f76e1feec7f59013448",
- "sha256:1b0ece94018ae21163d1f651b527156e1f03943b986188dd81bc7e066eae9d1c",
- "sha256:2e40be731ad618cb4974d5ba60d373cdf4f1b8dcbf1dcf4d9dff5e212baf69c5",
- "sha256:4ba59db1fcc27ea31368af524dcf874d9277f21fd2e1f7f1e2e0c75ee61419ed",
- "sha256:59ca9c6592da581a03d42cc4e270732552243dc45e87248aa8d636d53812f6a5",
- "sha256:5e0feb76849ca3e83dd396254e47c7dba65b3fa9ed3df67c2556293ae3e16de3",
- "sha256:6d205249a0293e62bbb3898c4c2e1ff8a22f98375a34775a259a0523111a8f6c",
- "sha256:6fcc5a3990e269f86d388f165a089259893851437b904f422d301cdce4ff25c8",
- "sha256:82847f2765835c8e5308f136bc34018d09b49037ec23ecc42b246424c767056b",
- "sha256:87902e5c03355335fc5992a74ba0247a70d937f326d852fc613b7f53516c0963",
- "sha256:9ab21d1cb156a620d3999dd92f7d1c86824c622873841d6b080ca5495fa10fef",
- "sha256:a1baa1dc8ecd88fb2d2a651671a84b9938461e8a8eed13e2f0a812a94084d1fa",
- "sha256:a244f7af80dacf21054386539699ce29bcc64796ed9850c99a34b41305630286",
- "sha256:a35af656a7ba1d3decdd4fae5322b87277de8ac98b7d9da657d9e212ece76a61",
- "sha256:b1fe1a6f3a6f355f6c29789b5927f8bd4f134a4bd9a781099a7c4f66af8850f5",
- "sha256:b5ad0adb51b2dee7d0ee75a69e9871e2ddfb061c73ea8bc439376298141f77f5",
- "sha256:ba3c7a2814ec8a176bb71f91478293d633c08582119e713a0c5351c0f77698da",
- "sha256:cd77d58fb2acf57c1d1ee2835567cd70e6f1835e32090538f17f8a3a99e5e34b",
- "sha256:cdb3a70285e8220875e4d2bc394e49b4988bdb1298ffa4e0bd81b2f613be397c",
- "sha256:deb529c40c3f1e38d53d5ae6cd077c21f1d49e13afc7936f7f868455e16b64a0",
- "sha256:e7894793e6e8540dbeac77c87b489e331947813511108ae097f1715c018b8f3d"
+ "sha256:0aa2b318cf81eb1693fcfcbb8007e95e231d7e1aa24288137f3b19905736c3ee",
+ "sha256:163c78c04f47f26ca1b21068cea25ed7c5ecafe5f5ab2ea4895656a750582b56",
+ "sha256:1e37626bcb8895c4b3873fcfd54e9bfc5ffec8d0f525651d6985fcc5c6b6003c",
+ "sha256:264fd15590b3f02a1fbc095e7e1f37cdac698ff3829e12ffdcffdce3772f9d44",
+ "sha256:3d9e1554cd9b5999070c467b18e5ae3ebd7369f02706a8850816f576a954295f",
+ "sha256:40c24960cd5cec55222963f255858a1c47c6fa50a65a5b03fd7de75e3700eaaa",
+ "sha256:46f404314dbec78cb342904f9596f25f9b16e7cf304030f1339e553c8e77f51c",
+ "sha256:4847f0c993298b82fad809ea2916d857d0073dc17b0510fbbced663b3265929d",
+ "sha256:48e15612a8357393d176638c8f68a19273676877caea983f8baf188bad430379",
+ "sha256:6725d2797c65598778409aba8cd67077bb089d5b7d3d87c2719b206dc84ec05e",
+ "sha256:99f0ba97e369f02a21bb95faa3a0de55991fd5f0ece2e30a9e2eaebeac238921",
+ "sha256:a41f303b3f9157a31ce7203e3ca757a0c40c96669e72d9b6ee1bce8507638970",
+ "sha256:a4305564e93f5c4584f6758149fd446df39fd1e0a8c89ca0deb3cce56106a027",
+ "sha256:a551d8cc267c634774830086da42e4ba157fa41dd3b93982bc9501b284b0c689",
+ "sha256:a6bc9432c2640b008d5f29bad737714eb3e14bb8854878eacf3d7955c4e91c36",
+ "sha256:c60175d011a2e551a2f74c84e21e7c982489b96b6a5e4b030ecdeacf2914da68",
+ "sha256:e46e2384209c91996d5ec16744234d1c906ab79a701ce1a26155c9ec890b8dc8",
+ "sha256:e607b8cdc2ae5d5a63cd1bec30a15b5ed583ac6a39f04b7ba0f03fcfbf29c05b",
+ "sha256:e94a39d5c40fffe7696009dbd11bc14a349b377e03a384ed011e03d698787dd3",
+ "sha256:eb2286249ebfe8fcb5b425e5ec77e4736d53ee56d3ad296f8947f67150f495e3",
+ "sha256:fdee7540d12519865b423af411bd60ddb513d2eb2cd921149b732854995bbf8b"
],
"index": "pypi",
- "version": "==1.18.2"
+ "version": "==1.18.3"
},
"obspy": {
"hashes": [
@@ -562,20 +562,38 @@
],
"version": "==0.13.2"
},
+ "typing": {
+ "hashes": [
+ "sha256:91dfe6f3f706ee8cc32d38edbbf304e9b7583fb37108fef38229617f8b3eba23",
+ "sha256:c8cabb5ab8945cd2f54917be357d134db9cc1eb039e59d1606dc1e60cb1d9d36",
+ "sha256:f38d83c5a7a7086543a0f649564d661859c5146a85775ab90c0d2f93ffaa9714"
+ ],
+ "index": "pypi",
+ "version": "==3.7.4.1"
+ },
+ "typing-extensions": {
+ "hashes": [
+ "sha256:6e95524d8a547a91e08f404ae485bbb71962de46967e1b71a0cb89af24e761c5",
+ "sha256:79ee589a3caca649a9bfd2a8de4709837400dfa00b6cc81962a1e6a1815969ae",
+ "sha256:f8d2bd89d25bc39dabe7d23df520442fa1d8969b82544370e03d88b5a591c392"
+ ],
+ "index": "pypi",
+ "version": "==3.7.4.2"
+ },
"urllib3": {
"hashes": [
- "sha256:2f3db8b19923a873b3e5256dc9c2dedfa883e33d87c690d9c7913e1f40673cdc",
- "sha256:87716c2d2a7121198ebcb7ce7cccf6ce5e9ba539041cfbaeecfb641dc0bf6acc"
+ "sha256:3018294ebefce6572a474f0604c2021e33b3fd8006ecd11d62107a5d2a963527",
+ "sha256:88206b0eb87e6d677d424843ac5209e3fb9d0190d0ee169599165ec25e9d9115"
],
- "version": "==1.25.8"
+ "version": "==1.25.9"
},
"uvicorn": {
"hashes": [
- "sha256:0f58170165c4495f563d8224b2f415a0829af0412baa034d6f777904613087fd",
- "sha256:6fdaf8e53bf1b2ddf0fe9ed06079b5348d7d1d87b3365fe2549e6de0d49e631c"
+ "sha256:a19de5b2f8dec56ec95bbbf3ce334b168efca6a238b2ba7c6cb020e3a42ce71a",
+ "sha256:aa6db2200cff18f8c550c869fab01f78ca987464d499f654e7aeb166fe0c0616"
],
"index": "pypi",
- "version": "==0.11.3"
+ "version": "==0.11.4"
},
"uvloop": {
"hashes": [
@@ -660,10 +678,10 @@
},
"click": {
"hashes": [
- "sha256:8a18b4ea89d8820c5d0c7da8a64b2c324b4dabb695804dbfea19b9be9d88c0cc",
- "sha256:e345d143d80bf5ee7534056164e5e112ea5e22716bbb1ce727941f4c8b471b9a"
+ "sha256:d2b5255c7c6349bc1bd1e59e08cd12acbbd63ce649f2588755783aa94dfb6b1a",
+ "sha256:dacca89f4bfadd5de3d7489b7c8a566eee0d3676333fbb50030263894c38c0dc"
],
- "version": "==7.1.1"
+ "version": "==7.1.2"
},
"coverage": {
"hashes": [
@@ -716,10 +734,10 @@
},
"identify": {
"hashes": [
- "sha256:2bb8760d97d8df4408f4e805883dad26a2d076f04be92a10a3e43f09c6060742",
- "sha256:faffea0fd8ec86bb146ac538ac350ed0c73908326426d387eded0bcc9d077522"
+ "sha256:23c18d97bb50e05be1a54917ee45cc61d57cb96aedc06aabb2b02331edf0dbf0",
+ "sha256:88ed90632023e52a6495749c6732e61e08ec9f4f04e95484a5c37b9caf40283c"
],
- "version": "==1.4.14"
+ "version": "==1.4.15"
},
"importlib-metadata": {
"hashes": [
@@ -765,11 +783,11 @@
},
"pre-commit": {
"hashes": [
- "sha256:487c675916e6f99d355ec5595ad77b325689d423ef4839db1ed2f02f639c9522",
- "sha256:c0aa11bce04a7b46c5544723aedf4e81a4d5f64ad1205a30a9ea12d5e81969e1"
+ "sha256:979b53dab1af35063a483bfe13b0fcbbf1a2cf8c46b60e0a9a8d08e8269647a1",
+ "sha256:f3e85e68c6d1cbe7828d3471896f1b192cfcf1c4d83bf26e26beeb5941855257"
],
"index": "pypi",
- "version": "==2.2.0"
+ "version": "==2.3.0"
},
"py": {
"hashes": [
@@ -892,10 +910,10 @@
},
"virtualenv": {
"hashes": [
- "sha256:00cfe8605fb97f5a59d52baab78e6070e72c12ca64f51151695407cc0eb8a431",
- "sha256:c8364ec469084046c779c9a11ae6340094e8a0bf1d844330fc55c1cefe67c172"
+ "sha256:5021396e8f03d0d002a770da90e31e61159684db2859d0ba4850fbea752aa675",
+ "sha256:ac53ade75ca189bc97b6c1d9ec0f1a50efe33cbf178ae09452dcd9fd309013c1"
],
- "version": "==20.0.17"
+ "version": "==20.0.18"
},
"waitress": {
"hashes": [
@@ -920,11 +938,11 @@
},
"webtest": {
"hashes": [
- "sha256:71114cd778a7d7b237ec5c8a5c32084f447d869ae62e48bcd5b73af211133e74",
- "sha256:da9cf14c103ff51a40dee4cac7657840d1317456eb8f0ca81289b5cbff175f4b"
+ "sha256:44ddfe99b5eca4cf07675e7222c81dd624d22f9a26035d2b93dc8862dc1153c6",
+ "sha256:aac168b5b2b4f200af4e35867cf316712210e3d5db81c1cbdff38722647bb087"
],
"index": "pypi",
- "version": "==2.0.34"
+ "version": "==2.0.35"
},
"zipp": {
"hashes": [
@@ -934,5 +952,4 @@
"version": "==3.1.0"
}
}
- }
}
diff --git a/geomagio/residual/Calculation.py b/geomagio/residual/Calculation.py
index e1494043c..8e2e2826d 100644
--- a/geomagio/residual/Calculation.py
+++ b/geomagio/residual/Calculation.py
@@ -1,68 +1,198 @@
from typing import List, Tuple
+from typing_extensions import Literal
import numpy as np
from pydantic import BaseModel
from .Absolute import Absolute
from .Angle import from_dms, to_dms
-from .MeasurementType import MeasurementType as mt
+from .MeasurementType import (
+ MeasurementType as mt,
+ DECLINATION_TYPES,
+ INCLINATION_TYPES,
+ MARK_TYPES,
+)
from .Measurement import (
AverageMeasurement,
Measurement,
average_measurement,
measurement_index,
)
+from .Reading import Reading
-# specify mark measurement types
-MARK_TYPES = [
- mt.FIRST_MARK_DOWN,
- mt.FIRST_MARK_UP,
- mt.SECOND_MARK_DOWN,
- mt.SECOND_MARK_UP,
-]
-
-# define measurement types used to calculate inclination
-INCLINATION_TYPES = [mt.NORTH_DOWN, mt.NORTH_UP, mt.SOUTH_DOWN, mt.SOUTH_UP]
+def calculate(reading: Reading) -> Reading:
+ """Calculate absolutes and scale value using residual method.
-# define measurement types used to calculate declination
-DECLINATION_TYPES = [mt.EAST_UP, mt.EAST_DOWN, mt.WEST_UP, mt.WEST_DOWN]
+ Parameters
+ -------
+ reading: reading to calculate absolutes from.
-
-def calculate(reading) -> Tuple[Absolute, Absolute, Absolute]:
- """
- 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.
+ Returns
+ -------
+ new reading object with calculated absolutes and scale_value.
+ NOTE: rest of reading object is shallow copy.
"""
- # gather oridinates, measuremenets, and metadata objects from reading
- index = measurement_index(reading.measurements)
# reference measurement, used to adjust absolutes
- reference = index[mt.WEST_DOWN][0]
+ reference = reading[mt.WEST_DOWN][0]
# calculate inclination
- inclination, f, mean = calculate_I(reading)
+ inclination, f, mean = calculate_I(
+ hemisphere=reading.hemisphere, measurements=reading.measurements
+ )
+ corrected_f = f + reading.pier_correction # TODO: should this be returned?
# calculate absolutes
- absoluteH, absoluteZ = calculate_HZ_absolutes(inclination, f, mean, reference)
- scale_value = calculate_scale(inclination, f, index[mt.NORTH_DOWN_SCALE])
- # calculate declination absolute and baseline
- absoluteD = calculate_D_absolute(reading, absoluteH, reference)
- return absoluteD, absoluteH, absoluteZ
+ absoluteH, absoluteZ = calculate_HZ_absolutes(
+ corrected_f=corrected_f, inclination=inclination, mean=mean, reference=reference
+ )
+ absoluteD = calculate_D_absolute(
+ azimuth=reading.azimuth,
+ h_baseline=absoluteH.baseline,
+ measurements=reading.measurements,
+ reference=reference,
+ )
+ # calculate scale
+ scale_value = calculate_scale_value(
+ corrected_f=corrected_f,
+ inclination=inclination,
+ measurements=reading[mt.NORTH_DOWN_SCALE],
+ )
+ # create new reading object
+ calculated = Reading(
+ absolutes=[absoluteD, absoluteH, absoluteZ],
+ scale_value=scale_value,
+ # copy other attributes
+ **reading.dict(exclude={"absolutes", "scale_value"}),
+ )
+ return calculated
+
+def calculate_D_absolute(
+ measurements: List[Measurement],
+ azimuth: float,
+ h_baseline: float,
+ reference: Measurement,
+) -> Absolute:
+ """Calculate D absolute.
-def calculate_I(reading) -> Tuple[float, float, Measurement]:
+ Parameters
+ ----------
+ measurements: list with at least declination and mark measurements.
+ azimuth: azimuth to mark in decimal degrees.
+ h_baseline: calculated H baseline value.
+ reference: reference measurement used to adjust absolute.
+
+ Returns
+ -------
+ D Absolute
"""
- Calculate inclination angles from measurements, ordinates,
- average ordinates from every measurement, and metadata.
- Returns inclination angle and calculated average f
+ # mean across all declination measurements
+ mean = average_measurement(measurements, DECLINATION_TYPES)
+ # average mark
+ average_mark = average_measurement(measurements, MARK_TYPES).angle
+ # adjust based on which is larger
+ mark_up = average_measurement(measurements, [mt.FIRST_MARK_UP]).angle
+ mark_down = average_measurement(measurements, [mt.FIRST_MARK_DOWN]).angle
+ if mark_up < mark_down:
+ average_mark += 90
+ else:
+ average_mark -= 90
+ # declination measurements
+ declination_measurements = [
+ average_measurement(measurements, [t]) for t in DECLINATION_TYPES
+ ]
+ # average declination meridian
+ meridian = np.average(
+ [
+ m.angle
+ + np.degrees(
+ m.measurement_type.meridian
+ * (np.arcsin(m.residual / np.sqrt((m.h + h_baseline) ** 2 + m.e ** 2)))
+ )
+ - np.degrees(np.arctan(m.e / (m.h + h_baseline)))
+ for m in declination_measurements
+ ]
+ )
+ # add subtract average mark angle from average meridian angle and add
+ # azimuth to get the declination baseline
+ d_b = (meridian - average_mark) + azimuth
+ # calculate absolute
+ d_abs = d_b + np.degrees(np.arctan(reference.e / (reference.h + h_baseline)))
+ return Absolute(element="D", absolute=d_abs, baseline=d_b)
+
+
+def calculate_HZ_absolutes(
+ inclination: float,
+ corrected_f: float,
+ mean: AverageMeasurement,
+ reference: Measurement,
+) -> Tuple[Absolute, Absolute]:
+ """Calculate H and Z absolutes.
+
+ Parameters
+ ----------
+ inclination: calculated inclination.
+ corrected_f: calculated f with pier correction.
+ mean: mean of inclination ordinates.
+ reference: reference measurement used to adjust absolutes.
+
+ Returns
+ -------
+ Tuple
+ - H Absolute
+ - Z Absolute
+ """
+ inclination_radians = np.radians(inclination)
+ h_abs = corrected_f * np.cos(inclination_radians)
+ z_abs = corrected_f * np.sin(inclination_radians)
+ h_b = round(np.sqrt(h_abs ** 2 - mean.e ** 2) - mean.h, 1)
+ z_b = round(z_abs - mean.z, 1)
+ # adjust absolutes to reference measurement
+ h_abs = np.sqrt((h_b + reference.h) ** 2 + (reference.e) ** 2)
+ z_abs = z_b + reference.z
+ # return absolutes
+ return (
+ Absolute(
+ element="H",
+ baseline=h_b,
+ absolute=h_abs,
+ starttime=mean.time,
+ endtime=mean.endtime,
+ ),
+ Absolute(
+ element="Z",
+ baseline=z_b,
+ absolute=z_abs,
+ starttime=mean.time,
+ endtime=mean.endtime,
+ ),
+ )
+
+
+def calculate_I(
+ measurements: List[Measurement], hemisphere: Literal[-1, 1] = 1
+) -> Tuple[float, float, Measurement]:
+ """Calculate inclination and f from measurements.
+
+ Parameters
+ ----------
+ measurements: list with at least inclination measurements.
+ hemisphere: +1 for northern hemisphere (default), -1 for southern hemisphere.
+
+ Returns
+ -------
+ Tuple
+ - inclination angle in decimal degrees
+ - uncorrected calculated f
+ - mean of inclination measurements
"""
- index = measurement_index(reading.measurements)
# mean across all inclination measurements
- mean = reading.get_average(INCLINATION_TYPES)
+ mean = average_measurement(measurements, INCLINATION_TYPES)
# mean within each type
- measurements = [reading.get_average(t) for t in INCLINATION_TYPES]
+ inclination_measurements = [
+ average_measurement(measurements, [t]) for t in INCLINATION_TYPES
+ ]
# get initial inclination angle, assumed to be the southdown angle
- inclination = reading.get_average([mt.NORTH_UP]).angle
+ inclination = average_measurement(measurements, [mt.NORTH_DOWN]).angle
if inclination >= 90:
inclination -= 180
# loop until inclination converges
@@ -72,7 +202,7 @@ def calculate_I(reading) -> Tuple[float, float, Measurement]:
last_inclination = inclination
# Update measurement f based on inclination
inclination_radians = np.radians(inclination)
- for m in measurements:
+ for m in inclination_measurements:
m.f = (
mean.f
+ (m.h - mean.h) * np.cos(inclination_radians)
@@ -88,111 +218,45 @@ def calculate_I(reading) -> Tuple[float, float, Measurement]:
* (
+m.angle
+ m.measurement_type.direction
- * (reading.hemisphere * np.degrees(np.arcsin(m.residual / m.f)))
+ * (hemisphere * np.degrees(np.arcsin(m.residual / m.f)))
)
)
- for m in measurements
+ for m in inclination_measurements
]
)
- # calculate corrected f
- f = np.average([m.f for m in measurements]) + reading.pier_correction
+ # calculate uncorrected f
+ f = np.average([m.f for m in inclination_measurements])
return inclination, f, mean
-def calculate_HZ_absolutes(
- inclination: float, f: float, mean: AverageMeasurement, reference: Measurement
-):
- inclination_radians = np.radians(inclination)
- h_abs = f * np.cos(inclination_radians)
- z_abs = f * np.sin(inclination_radians)
- h_b = round(np.sqrt(h_abs ** 2 - mean.e ** 2) - mean.h, 1)
- z_b = round(z_abs - mean.z, 1)
- # adjust absolutes to reference measurement
- h_abs = np.sqrt((h_b + reference.h) ** 2 + (reference.e) ** 2)
- z_abs = z_b + reference.z
- # return absolutes
- return (
- Absolute(
- element="H",
- baseline=h_b,
- absolute=h_abs,
- starttime=mean.time,
- endtime=mean.endtime,
- ),
- Absolute(
- element="Z",
- baseline=z_b,
- absolute=z_abs,
- starttime=mean.time,
- endtime=mean.endtime,
- ),
- )
-
-
-def calculate_D_absolute(
- reading, absoluteH: Absolute, reference: Measurement
-) -> Absolute:
- """
- Calculate declination absolute and declination baseline from
- ordinates, measurements, measurement_index(dictionary), azimuth and H baseline
- Returns absolute and baseline for declination.
- """
- index = measurement_index(reading.measurements)
- # mean across all declination measurements
- mean = reading.get_average(DECLINATION_TYPES)
- # average mark
- average_mark = reading.get_average(MARK_TYPES).angle
- # adjust based on which is larger
- mark_up = index[mt.FIRST_MARK_UP][0]
- mark_down = index[mt.FIRST_MARK_DOWN][0]
- if mark_up.angle < mark_down.angle:
- average_mark += 90
- else:
- average_mark -= 90
- # declination measurements
- measurements = [reading.get_average(t) for t in DECLINATION_TYPES]
- # average declination meridian
- h_base = absoluteH.baseline
- meridian = np.average(
- [
- (
- m.angle
- + np.degrees(
- m.measurement_type.meridian
- * (np.arcsin(m.residual / np.sqrt((m.h + h_base) ** 2 + m.e ** 2)))
- )
- - np.degrees(np.arctan(m.e / (m.h + h_base)))
- )
- for m in measurements
- ]
- )
- # add subtract average mark angle from average meridian angle and add
- # azimuth to get the declination baseline
- d_b = (meridian - average_mark) + reading.azimuth
- # calculate absolute
- d_abs = d_b + np.degrees(np.arctan(reference.e / (reference.h + h_base)))
- return Absolute(element="D", absolute=d_abs, baseline=d_b)
+def calculate_scale_value(
+ measurements: List[Measurement], inclination: float, corrected_f: float
+) -> float:
+ """Calculate scale value.
+ Parameters
+ ----------
+ measurements: measurements to be used for scale calculation.
+ should have type NORTH_DOWN_SCALE.
+ inclination: calculated inclination.
+ corrected_f: calculated f with pier correction.
-def calculate_scale(inclination: float, f: float, measurements: List[Measurement]):
- """
- Calculate scale value from calulated f(from inclination computations),
- 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.
+ Returns
+ -------
+ Calculated scale value.
"""
inclination_radians = np.radians(inclination)
- scale0, scale1 = measurements[0], measurements[-1]
+ m1, m2 = measurements[0], measurements[-1]
field_change = (
np.degrees(
(
- -np.sin(inclination_radians) * (scale1.h - scale0.h)
- + np.cos(inclination_radians) * (scale1.z - scale0.z)
+ -np.sin(inclination_radians) * (m2.h - m1.h)
+ + np.cos(inclination_radians) * (m2.z - m1.z)
)
- / f
+ / corrected_f
)
+ 0.1668
)
- residual_change = scale1.residual - scale0.residual
- scale_value = f * field_change / np.abs(residual_change)
+ residual_change = m2.residual - m1.residual
+ scale_value = corrected_f * field_change / np.abs(residual_change)
return scale_value
diff --git a/geomagio/residual/Measurement.py b/geomagio/residual/Measurement.py
index 1bad924cf..86fcd7403 100644
--- a/geomagio/residual/Measurement.py
+++ b/geomagio/residual/Measurement.py
@@ -34,15 +34,24 @@ class AverageMeasurement(Measurement):
endtime: Optional[UTCDateTime] = None
-def average_measurement(measurements: List[Measurement]) -> AverageMeasurement:
+def average_measurement(
+ measurements: List[Measurement], types: List[MeasurementType] = None
+) -> AverageMeasurement:
"""Calculate average from multiple measurements.
- returns None if measurements is empty or None
- otherwise returns Measurement
- - using type from first measurement,
- - with empty time,
- - averaging all other values
+ Parameters
+ ----------
+ measurements - source measurements for average
+ types - optional list of types to include, default all
+
+ Returns
+ -------
+ None - if no measurements
+ Otherwise, average of matching measurements.
+ Type is copied from first measurement.
"""
+ if types:
+ measurements = [m for m in measurements if m.measurement_type in types]
if len(measurements) == 0:
# no measurements to average
return None
diff --git a/geomagio/residual/MeasurementType.py b/geomagio/residual/MeasurementType.py
index 2b43b3b14..f8e1216ae 100644
--- a/geomagio/residual/MeasurementType.py
+++ b/geomagio/residual/MeasurementType.py
@@ -1,6 +1,6 @@
import enum
-# FIXME: add another measurement type that is a scaling measurement type
+
class MeasurementType(str, enum.Enum):
"""Measurement types used during absolutes."""
@@ -35,17 +35,6 @@ class MeasurementType(str, enum.Enum):
else:
return -1
- @property
- def shift(self):
- if self == MeasurementType.SOUTH_DOWN:
- return -180
- if self == MeasurementType.SOUTH_UP:
- return 180
- if self == MeasurementType.NORTH_UP:
- return 0
- if self == MeasurementType.NORTH_DOWN:
- return 360
-
@property
def meridian(self):
if self in [
@@ -57,3 +46,41 @@ class MeasurementType(str, enum.Enum):
return 1
else:
return -1
+
+ @property
+ def shift(self):
+ if self == MeasurementType.SOUTH_DOWN:
+ return -180
+ if self == MeasurementType.SOUTH_UP:
+ return 180
+ if self == MeasurementType.NORTH_UP:
+ return 0
+ if self == MeasurementType.NORTH_DOWN:
+ return 360
+
+
+# define measurement types used to calculate declination
+DECLINATION_TYPES = [
+ MeasurementType.EAST_UP,
+ MeasurementType.EAST_DOWN,
+ MeasurementType.WEST_UP,
+ MeasurementType.WEST_DOWN,
+]
+
+
+# define measurement types used to calculate inclination
+INCLINATION_TYPES = [
+ MeasurementType.NORTH_DOWN,
+ MeasurementType.NORTH_UP,
+ MeasurementType.SOUTH_DOWN,
+ MeasurementType.SOUTH_UP,
+]
+
+
+# specify mark measurement types
+MARK_TYPES = [
+ MeasurementType.FIRST_MARK_DOWN,
+ MeasurementType.FIRST_MARK_UP,
+ MeasurementType.SECOND_MARK_DOWN,
+ MeasurementType.SECOND_MARK_UP,
+]
diff --git a/geomagio/residual/Reading.py b/geomagio/residual/Reading.py
index be6b90c70..19041231f 100644
--- a/geomagio/residual/Reading.py
+++ b/geomagio/residual/Reading.py
@@ -1,12 +1,12 @@
import collections
from typing import Dict, List, Optional
+from typing_extensions import Literal
from pydantic import BaseModel
from .Absolute import Absolute
from .Measurement import AverageMeasurement, Measurement, average_measurement
from .MeasurementType import MeasurementType
-from .Calculation import calculate
class Reading(BaseModel):
@@ -22,18 +22,17 @@ class Reading(BaseModel):
pier_correction: pier correction value, nT.
"""
- absolutes: Optional[List[Absolute]] = None
+ absolutes: List[Absolute] = []
azimuth: float = 0
- hemisphere: float = 1
- measurements: Optional[List[Measurement]] = []
- metadata: Optional[Dict] = []
+ hemisphere: Literal[-1, 1] = 1
+ measurements: List[Measurement] = []
+ metadata: Dict = {}
pier_correction: float = 0
+ scale_value: float = None
- def get_average(self, types: List[MeasurementType]) -> AverageMeasurement:
- return average_measurement(
- [m for m in self.measurements if m.measurement_type in types]
- )
+ def __getitem__(self, measurement_type: MeasurementType):
+ """Provide access to measurements by type.
- def update_absolutes(self):
- self.absolutes = calculate(self)
- return self.absolutes
+ Example: reading[MeasurementType.WEST_DOWN]
+ """
+ return [m for m in self.measurements if m.measurement_type == measurement_type]
diff --git a/geomagio/residual/__init__.py b/geomagio/residual/__init__.py
index 9ec21abb7..59e03297a 100644
--- a/geomagio/residual/__init__.py
+++ b/geomagio/residual/__init__.py
@@ -1,11 +1,23 @@
# residual module
from __future__ import absolute_import
-from .Absolute import Absolute
from . import Angle
+from .Absolute import Absolute
+from .Calculation import (
+ calculate,
+ calculate_D_absolute,
+ calculate_HZ_absolutes,
+ calculate_I,
+ calculate_scale_value,
+)
from .CalFileFactory import CalFileFactory
-from .Measurement import Measurement
-from .MeasurementType import MeasurementType
+from .Measurement import Measurement, AverageMeasurement, average_measurement
+from .MeasurementType import (
+ MeasurementType,
+ DECLINATION_TYPES,
+ INCLINATION_TYPES,
+ MARK_TYPES,
+)
from .Reading import Reading
from .SpreadsheetAbsolutesFactory import SpreadsheetAbsolutesFactory
from .WebAbsolutesFactory import WebAbsolutesFactory
@@ -13,7 +25,16 @@ from .WebAbsolutesFactory import WebAbsolutesFactory
__all__ = [
"Absolute",
"Angle",
- "CalFileFactory",
+ "AverageMeasurement",
+ "average_measurement" "CalFileFactory",
+ "calculate",
+ "calculate_D_absolute",
+ "calculate_HZ_absolutes",
+ "calculate_I",
+ "calculate_scale_value",
+ "DECLINATION_TYPES",
+ "INCLINATION_TYPES",
+ "MARK_TYPES",
"Measurement",
"MeasurementType",
"Reading",
diff --git a/test/residual_test/residual_test.py b/test/residual_test/residual_test.py
index 0266ad318..c132d48ed 100644
--- a/test/residual_test/residual_test.py
+++ b/test/residual_test/residual_test.py
@@ -1,7 +1,8 @@
-from geomagio.residual import SpreadsheetAbsolutesFactory
from numpy.testing import assert_almost_equal
import pytest
+from geomagio.residual import calculate, SpreadsheetAbsolutesFactory
+
def assert_absolutes(original, result):
"""
@@ -37,8 +38,8 @@ def compare_spreadsheet_absolutes(path):
reading = saf.parse_spreadsheet(path=path)
# establish original absolute object
original = {a.element: a for a in reading.absolutes}
- reading.update_absolutes()
- result = {a.element: a for a in reading.absolutes}
+ calculated = calculate(reading)
+ result = {a.element: a for a in calculated.absolutes}
return original, result
--
GitLab