Mandatory REL file to draw SRF planes, check if fully inside VM domain (

#238)

* Optional REL file to draw SRF planes, check if fully inside VM domain

* Added a comment

* Use NZTM coords and validate_vm_bounds()

* Removed unintended changes

* Changes accomodating feedback

VM/plot_vm.py

@@ -1,12 +1,17 @@
+"""
+Plots the VM domain and the SRF planes
+"""
+
 from argparse import ArgumentParser
+from io import StringIO
 from logging import Logger
 from pathlib import Path
 from tempfile import TemporaryDirectory
 
 import numpy as np
 import yaml
 
-from qcore import geo, gmt, qclogging
+from qcore import geo, gmt, qclogging, validate_vm
 
 
 def plot_vm(
@@ -24,23 +29,37 @@ def plot_vm(
 
     Parameters
     ----------
-    vm_params_dict :
-    srf_corners :
-    land_outline_path :
-    centre_line_path :
-    mag :
-    outdir :
-    ptemp :
-    logger :
+    vm_params_dict : dict
+        Dictionary extracted from vm_params.yaml file.
+    srf_corners : np.ndarray.
+        Corners of SRF planes. Formatted as [[[lon,lat],[lon,lat],[lon,lat],[lon,lat]],...]
+    land_outline_path : Path
+        Path to the land outline file
+    centre_line_path : Path
+        Path to the centre line file
+    mag : float
+        Magnitude of the event
+    outdir : Path
+        Output directory
+    ptemp : Path
+        Temporary directory
+    logger : Logger
+        Default is qclogging.get_basic_logger()
     """
 
+    from rel2vm_params import write_srf_path
+
     logger.debug("Plotting vm")
     p = gmt.GMTPlot(ptemp / "optimisation.ps")
     p.spacial("M", vm_params_dict["plot_region"], sizing=7)
     p.coastlines()
 
-    srf_path = ptemp / "srf.path"
-    if srf_path.exists():
+    # SRF domain
+    if len(srf_corners) > 0:
+        srf_path = write_srf_path(
+            srf_corners, ptemp
+        )  # write srf.path file if not already present, otherwise use the existing file
+
         # filled slip area
         p.path(srf_path, is_file=True, fill="yellow", split="-")
         # top edge
@@ -49,7 +68,7 @@ def plot_vm(
                 "\n".join([" ".join(map(str, ll)) for ll in plane[:2]]), is_file=False
             )
 
-    # vm domain (simple and adjusted)
+    # plot the VM domain (simple and adjusted)
     p.path(vm_params_dict["path"], is_file=False, close=True, fill="black@95")
     if vm_params_dict["adjusted"]:
         p.path(
@@ -119,24 +138,31 @@ def main(
     name: str,
     vm_params_dict: dict,
     outdir: Path,
+    rel_path: Path,
     logger: Logger = qclogging.get_basic_logger(),
 ):
     """
-    vm_params_dict loaded from vm_params.yaml doesn't have all info plot_vm() needs.
-    This function gathers and works out the necessary input (except SRF-relevant info) to run this file as a stand-alone script
+    Gathers necessary input to call plot_vm() function to plot VM domain and SRF planes (if realisation CSV is supplied)
 
     Parameters
     ----------
-    name : name of the fault/event
-    vm_params_dict : Dictionary extracted from vm_params.yaml
-    outdir :
+    name : str
+        name of the fault/event. This is used to name the output file.
+    vm_params_dict : dict
+        Dictionary extracted from vm_params.yaml file.
+    outdir : Path
+        Output directory
+    rel_path : Path
+        Path to the realisation csv file.
+
     logger :
     """
     from rel2vm_params import (
         get_vm_land_proportion,
         corners2region,
         NZ_CENTRE_LINE,
         NZ_LAND_OUTLINE,
+        load_rel,
     )
 
     # vm_params_dict is the dictionary directly loaded from vm_params.yaml
@@ -173,10 +199,13 @@ def main(
 
         vm_params_dict["plot_region"] = plot_region
 
-        # plotting the domain of VM. No SRF
+        srf_meta = load_rel(rel_path)
+        srf_corners = srf_meta["corners"]
+
+        # plotting the domain of VM.
         plot_vm(
             vm_params_dict,
-            [],
+            srf_corners,
             NZ_LAND_OUTLINE,
             NZ_CENTRE_LINE,
             vm_params_dict["mag"],
@@ -185,25 +214,34 @@ def main(
             logger=logger,
         )
 
+        # Validate the VM domain. Only needed when this code is run as a standalone script.
+        # If this code is run as part of the VM workflow, the validation is done in the main script.
+        polygon = np.loadtxt(StringIO(vm_params_dict["path"]))
+        errors = validate_vm.validate_region(polygon)
+        errors.extend(validate_vm.validate_vm_bounds(polygon, srf_corners))
+        if errors:
+            logger.warning(f"WARNING: {errors}")
+
 
 def load_args(logger: Logger = qclogging.get_basic_logger()):
     """
     Unpacks arguments and does basic checks
 
     Parameters
     ----------
-    logger :
+    logger :  Logger
+        Default is qclogging.get_basic_logger()
 
     Returns
     -------
+    Processed arguments
 
     """
     parser = ArgumentParser()
     arg = parser.add_argument
 
-    arg("name", help="Name of the fault")
-
-    arg("vm_params_path", help="path to vm_params.yaml")
+    arg("vm_params_path", help="path to vm_params.yaml", type=Path)
+    arg("rel_path", help="Path to the realisation csv file", type=Path)
 
     arg(
         "-o",
@@ -214,12 +252,22 @@ def load_args(logger: Logger = qclogging.get_basic_logger()):
     )
 
     args = parser.parse_args()
-    args.vm_params_path = Path(args.vm_params_path).resolve()
+    args.vm_params_path = args.vm_params_path.resolve()
+    args.rel_path = args.rel_path.resolve()
 
-    if args.outdir is None:
+    if (
+        args.outdir is None
+    ):  # if not specified, use the directory that contains vm_params.yaml
         args.outdir = args.vm_params_path.parent
+
     args.outdir = Path(args.outdir).resolve()
 
+    args.outdir.mkdir(exist_ok=True, parents=True)
+    assert args.vm_params_path.exists(), f"File is not present: {args.vm_params_path}"
+    assert args.rel_path.exists(), f"File is not present: {args.rel_path}"
+
+    args.name = args.rel_path.stem
+
     return args
 
 
@@ -231,4 +279,4 @@ def load_args(logger: Logger = qclogging.get_basic_logger()):
     with open(args.vm_params_path, "r") as f:
         vm_params_dict = yaml.load(f, Loader=yaml.SafeLoader)
 
-    main(args.name, vm_params_dict, args.outdir, logger=logger)
+    main(args.name, vm_params_dict, args.outdir, args.rel_path, logger=logger)

VM/rel2vm_params.py

@@ -22,19 +22,18 @@
 import sys
 from tempfile import TemporaryDirectory
 
-from srf_generation.input_file_generation.realisation_to_srf import get_corners_dbottom
 
 from qcore import constants, geo, gmt, qclogging
 from qcore.geo import R_EARTH
 from qcore.utils import dump_yaml
 from qcore.simulation_structure import get_fault_from_realisation
-from qcore.validate_vm import validate_vm_bounds
-
+from qcore.validate_vm import validate_vm_bounds, validate_region
+from srf_generation.input_file_generation.realisation_to_srf import get_corners_dbottom
 from VM.models.classdef import Site, Fault, TectType, estimate_z1p0, FaultStyle
 from VM.models.Bradley_2010_Sa import Bradley_2010_Sa
 from VM.models.AfshariStewart_2016_Ds import Afshari_Stewart_2016_Ds
 
-from plot_vm import plot_vm
+from VM.plot_vm import plot_vm
 
 script_dir = Path(__file__).resolve().parent
 NZ_CENTRE_LINE = script_dir / "../SrfGen/NHM/res/centre.txt"
@@ -559,6 +558,47 @@ def centre_lon(lat_target: float):
     return lon_target
 
 
+def write_srf_path(srf_corners: np.ndarray, wd: Path):
+    """
+    Write a temporary file srf.path containing SRF corner coordinates.
+    Used for plotting SRF planes using GMT. Needs to have the first point repeated at the end to close the path.
+
+    Example
+    -------
+    > srf plane
+    167.875309 -45.036523
+    167.498049 -45.811822
+    167.414717 -45.801548
+    167.793114 -45.026250
+    167.875309 -45.036523
+    > srf plane
+    167.496522 -45.811651
+    167.348478 -46.000020
+    167.264864 -45.989745
+    167.413191 -45.801377
+    167.496522 -45.811651
+
+
+    Parameters
+    ----------
+    srf_corners : corners of the srf
+    wd : working directory
+
+    Returns
+    -------
+    path to the srf path file
+
+    """
+    srf_path = wd / "srf.path"
+    if not srf_path.exists():
+        with open(srf_path, "w") as sp:
+            for plane in srf_corners:
+                sp.write("> srf plane\n")
+                np.savetxt(sp, plane, fmt="%.6f")
+                sp.write(f"{plane[0][0]:6f} {plane[0][1]:6f}\n")  # close the plane
+    return srf_path
+
+
 def optimise_vm_params(
     srf_meta: dict,
     ds_multiplier: float,
@@ -629,11 +669,7 @@ def optimise_vm_params(
     )
 
     # for plotting and calculating VM domain distance
-    with open(temp_dir / "srf.path", "wb") as sp:
-        for plane in srf_meta["corners"]:
-            sp.write("> srf plane\n".encode())
-            np.savetxt(sp, plane, fmt="%f")
-            sp.write("%f %f\n".encode() % (tuple(plane[0])))
+    write_srf_path(srf_meta["corners"], temp_dir)
 
     if fault_depth > rrup and not deep_rupture:
         logger.warning(
@@ -749,13 +785,19 @@ def optimise_vm_params(
         )
         success = False
 
-    bounds_invalid = validate_vm_bounds([c1, c2, c3, c4], srf_meta["corners"].tolist())
+    # check if corners are within NZ
+    bounds_invalid = validate_region([c1, c2, c3, c4])
+    bounds_invalid.extend(
+        validate_vm_bounds([c1, c2, c3, c4], srf_meta["corners"].tolist())
+    )
+
     if bounds_invalid:
         logger.warning(f"Bounds not valid, not making VM: {bounds_invalid}")
         success = False
 
     # modified sim time
     vm_corners = np.asarray([c1, c2, c3, c4])
+
     initial_time = get_sim_duration(
         vm_corners,
         np.concatenate(srf_meta["corners"], axis=0),