Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Typer CLI Interface #2

Merged
merged 9 commits into from
May 16, 2024
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,23 @@ This repository contains a script and schema to generate a database storing faul
- qcore
- sqlite3

## Setting up development environment

Use [Poetry](https://python-poetry.org/) to setup and manage the virtual enviroment. From within the repo run

```bash
$ poetry lock && poetry install
```

This will install all the dependencies required for development. To start a virtual environment run

```bash
$ poetry shell
```

To integrate with PyCharm (if you use it) see [PyCharm's Documentation](https://www.jetbrains.com/help/pycharm/poetry.html#poetry-pyproject)


## Generate your own database
1. Obtain the CRU fault system solution (ask a member of the software team for this file). It should be a zip file with the following structure
```
Expand Down
232 changes: 183 additions & 49 deletions nshm_geojson_fault_parser.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,38 +2,137 @@
import csv
import json
import sqlite3
from pathlib import Path
from sqlite3 import Connection
from typing import Any
from typing import Annotated, Any

import numpy as np
import qcore.geo
import typer

app = typer.Typer()

def strike_between_coordinates(a: (float, float), b: (float, float)) -> float:
a_lat, a_lon = a
b_lat, b_lon = b

def strike_between_coordinates(
lispandfound marked this conversation as resolved.
Show resolved Hide resolved
point_a_coords: (float, float), point_b_coords: (float, float)
) -> float:
"""Compute the strike angle between two points given as (lat, lon) pairs.

Parameters
----------
point_a_coords : (float, float)
The coordinates of the first point.
point_b_coords : (float, float)
The coordinates of the second point.

Returns
-------
float
The strike angle (in degrees) between point a and point b.

"""
a_lat, a_lon = point_a_coords
b_lat, b_lon = point_b_coords
return qcore.geo.ll_bearing(a_lon, a_lat, b_lon, b_lat)


def distance_between(a: (float, float), b: (float, float)) -> float:
a_lat, a_lon = a
b_lat, b_lon = b
def distance_between(
lispandfound marked this conversation as resolved.
Show resolved Hide resolved
point_a_coords: (float, float), point_b_coords: (float, float)
) -> float:
"""Compute the distance between two points given as (lat, lon) pairs.

Parameters
----------
point_a_coords : (float, float)
The coordinates of the first point.
point_b_coords : (float, float)
The coordinates of the second point.

Returns
-------
float
The distance (in kilometres) between point a and point b.

"""
a_lat, a_lon = point_a_coords
b_lat, b_lon = point_b_coords
return qcore.geo.ll_dist(a_lon, a_lat, b_lon, b_lat)


def centre_point(
a: (float, float), b: (float, float), dip: float, dip_dir: float, width: float
def centre_point_of_fault(
point_a_coords: (float, float),
point_b_coords: (float, float),
dip: float,
dip_dir: float,
width: float,
) -> (float, float):
a_lat, a_lon = a
b_lat, b_lon = b
"""Find the centre point of a fault defined by two points and a dip direction.

This function calculates the centre point of a fault given the coordinates
of two points on the fault trace, the dip angle, dip direction, and the
width of the fault. The centre point is defined as the midpoint of the
line connecting points a and a, shifted along the dip direction by half
the *projected* fault width. See the diagram.

point a point b
+───────────────+
│ │ │ │
│ │ │ │
│ │dip dir│ │
│ │ │ │
│ ∨ │ │ width
│ *centre │ │
│ │ │
│ │ │
│ │ │
└───────────────┘

The dip angle should be provided in degrees, with 0 degrees indicating
an (impossible) horizontal fault and 90 degrees indicating a vertical
fault. The dip direction is a bearing in degrees from north.

Parameters
----------
a : (float, float)
The (lat, lon) coordinates of the point a.
b : (float, float)
The (lat, lon) coordinates of the point b.
dip : float
The dip angle of the fault.
dip_dir : float
The dip direction of the fault.
width : float
The width of the fault. Note this is the actual width, not the
projected width.

Returns
-------
(float, float)
The (lat, lon) coordinates of the (projected) centre point of
the fault.
"""
a_lat, a_lon = point_a_coords
b_lat, b_lon = point_b_coords
c_lon, c_lat = qcore.geo.ll_mid(a_lon, a_lat, b_lon, b_lat)
projected_width = width * np.cos(np.radians(dip)) / 2
return qcore.geo.ll_shift(c_lat, c_lon, projected_width, dip_dir)
projected_width = width * np.cos(np.radians(dip))
return qcore.geo.ll_shift(c_lat, c_lon, projected_width / 2, dip_dir)


def insert_magnitude_frequency_distribution(
conn: Connection, magnitude_frequency_distribution: list[dict[str, float | str]]
):
"""Inserts magnitude-frequency distribution data into a database.

Parameters
----------
conn : sqlite3.Connection
A connection object to the SQLite database.
magnitude_frequency_distribution : list[dict[str, Union[float, str]]]
A list of dictionaries containing magnitude-frequency distribution
data. Each dictionary should have keys 'Section Index', representing
the fault section index, and keys representing magnitudes associated
with their respective annual occurrence rate.
"""
for section_distribution in magnitude_frequency_distribution:
segment_id = int(section_distribution["Section Index"])
for magnitude_key, probability_raw in section_distribution.items():
Expand All @@ -42,12 +141,27 @@ def insert_magnitude_frequency_distribution(
magnitude = float(magnitude_key)
probability = float(probability_raw)
conn.execute(
"INSERT INTO magnitude_frequency_distribution (fault_id, magnitude, probability) VALUES (?, ?, ?)",
"""
INSERT INTO magnitude_frequency_distribution (
fault_id, magnitude, probability
) VALUES (?, ?, ?)
""",
(segment_id, magnitude, probability),
)


def insert_faults(conn: Connection, fault_map: dict[str, Any]):
"""Inserts fault data into the database.

Parameters
----------
conn : sqlite3.Connection
A connection object to the SQLite database.
fault_map : list[dict[str, Any]]
A list of dictionaries containing fault data. Each dictionary should
represent a fault feature and contain necessary properties like dip,
rake, depth, etc.
"""
for feature in fault_map:
properties = feature["properties"]
dip = properties["DipDeg"]
Expand All @@ -62,7 +176,8 @@ def insert_faults(conn: Connection, fault_map: dict[str, Any]):
parent_id = properties["ParentID"]
parent_name = properties["ParentName"]
conn.execute(
"INSERT OR REPLACE INTO parent_fault (parent_id, name) VALUES (?, ?)",
"""INSERT OR REPLACE INTO parent_fault (parent_id, name)
VALUES (?, ?)""",
(parent_id, parent_name),
)
conn.execute(
Expand All @@ -72,11 +187,25 @@ def insert_faults(conn: Connection, fault_map: dict[str, Any]):
for i in range(len(leading_edge) - 1):
left = tuple(reversed(leading_edge[i]))
right = tuple(reversed(leading_edge[i + 1]))
c_lat, c_lon = centre_point(left, right, dip, dip_dir, dbottom)
c_lat, c_lon = centre_point_of_fault(left, right, dip, dip_dir, dbottom)
strike = strike_between_coordinates(left, right)
length = distance_between(left, right)
conn.execute(
"INSERT INTO fault_segment (strike, rake, dip, dtop, dbottom, length, width, dip_dir, clon, clat, fault_id) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)",
"""
INSERT INTO fault_segment (
strike,
rake,
dip,
dtop,
dbottom,
length,
width,
dip_dir,
clon,
clat,
fault_id
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
""",
(
strike,
rake,
Expand All @@ -93,32 +222,17 @@ def insert_faults(conn: Connection, fault_map: dict[str, Any]):
)


def subsection_parent_map(fault_features: dict[str, Any]) -> dict[str, Any]:
merged = {}
for feature in fault_features:
parent_name = feature["properties"]["ParentName"]
if parent_name not in merged:
merged[parent_name] = {
"properties": feature["properties"],
"geometry": feature["geometry"]["coordinates"],
}
else:
merged[parent_name]["geometry"].extend(
feature["geometry"]["coordinates"][1:]
)
return merged


def subsection_parent_lookup(fault_features: dict[str, Any]) -> dict[str, str]:
subsection_parent_lookup_table = {}
for feature in fault_features:
fault_id = feature["properties"]["FaultID"]
parent_id = feature["properties"]["ParentID"]
subsection_parent_lookup_table[fault_id] = parent_id
return subsection_parent_lookup


def insert_ruptures(conn: Connection, indices: dict[int, int]):
"""Inserts rupture data into the database.

Parameters
----------
conn : sqlite3.Connection
A connection object to the SQLite database.
indices : dict[int, int]
A dictionary containing rupture indices mapped to fault indices. If
indices[0] = 1, then the fault 1 is involved in rupture 0.
"""
for row in indices:
rupture_idx, fault_idx = [int(value) for value in row.values()]
conn.execute(
Expand All @@ -131,19 +245,39 @@ def insert_ruptures(conn: Connection, indices: dict[int, int]):
)


if __name__ == "__main__":
with open("fault_sections.geojson", "r", encoding="utf-8") as fault_file:
@app.command()
def main(
fault_sections_geojson_filepath: Annotated[
Path,
typer.Option(help="Fault sections geojson file", readable=True, exists=True),
] = "fault_sections.geojson",
fast_indices_filepath: Annotated[
Path, typer.Option(help="Fast indices csv file", readable=True, exists=True)
] = "fast_indices.csv",
mfds_filepath: Annotated[
Path, typer.Option(help="MFDS filepath", readable=True, exists=True)
] = "sub_seismo_on_fault_mfds.csv",
sqlite_db_path: Annotated[
Path, typer.Option(help="Output SQLite DB path", writable=True, exists=True)
] = "nshm2022.db",
):
"""Generate the NSHM2022 rupture data from a CRU system solution package."""
with open(fault_sections_geojson_filepath, "r", encoding="utf-8") as fault_file:
geojson_object = json.load(fault_file)
with open("fast_indices.csv", "r", encoding="utf-8") as csv_file_handle:
with open(fast_indices_filepath, "r", encoding="utf-8") as csv_file_handle:
csv_reader = csv.DictReader(csv_file_handle)
indices = list(csv_reader)
with open(
"sub_seismo_on_fault_mfds.csv", "r", encoding="utf-8"
) as mfds_file_handle:
with open(mfds_filepath, "r", encoding="utf-8") as mfds_file_handle:
csv_reader = csv.DictReader(mfds_file_handle)
magnitude_frequency_distribution = list(csv_reader)
with sqlite3.connect("nshm2022.db") as conn:
with sqlite3.connect(sqlite_db_path) as conn:
# To enforce foreign key constraints.
conn.execute("PRAGMA foreign_keys = 1")

insert_faults(conn, geojson_object["features"])
insert_ruptures(conn, indices)
insert_magnitude_frequency_distribution(conn, magnitude_frequency_distribution)


if __name__ == "__main__":
app()
2 changes: 2 additions & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ readme = "README.md"
python = "^3.12"
numpy = "*"
qcore = { git = "https://github.com/ucgmsim/qcore" }
typer = "*"

[tool.poetry.group.dev.dependencies]
python-lsp-server = "*"
Expand All @@ -19,6 +20,7 @@ rope = "*"
pyflakes = "*"
mccabe = "*"
pycodestyle = "*"
doq = "*"


[build-system]
Expand Down
Loading