Fix/thicknesses_

README.md

@@ -1,4 +1,4 @@
-# Map2Loop 3.0
+# Map2Loop 3.1
 
 Generate 3D geological model inputs from geological maps — a high-level implementation and extension of the original map2loop code developed by Prof. Mark Jessell at UWA. To see an example interactive model built with map2loop and LoopStructural, follow this link:
 

dependencies.txt

@@ -4,5 +4,5 @@ networkx
 owslib
 map2model
 hjson
-loopprojectfile==0.1.3
+loopprojectfile==0.1.4
 beartype

map2loop/_datasets/geodata_files/load_map2loop_data.py

@@ -0,0 +1,43 @@
+import geopandas
+import pkg_resources
+from osgeo import gdal
+
+def load_hamersley_geology():
+    """
+    Loads Hamersley geology data from a shapefile
+
+    Args:
+        path (str):
+            The path to the shapefile
+
+    Returns:
+        geopandas.GeoDataFrame: The geology data
+    """
+    stream = pkg_resources.resource_filename("map2loop", "/_datasets/geodata_files/hamersley/geology.geojson")
+    return geopandas.read_file(stream)
+
+
+def load_hamersley_structure():
+    """
+    Loads Hamersley structure data from a shapefile
+
+    Args:
+        path (str):
+            The path to the shapefile
+
+    Returns:
+        geopandas.GeoDataFrame: The structure data
+    """
+
+    path = pkg_resources.resource_filename("map2loop", "/_datasets/geodata_files/hamersley/structure.geojson")
+    return geopandas.read_file(path)
+
+def load_hamersley_dtm():
+    """
+    Load DTM data from a raster file
+
+    Returns:
+        gdal.Dataset: The DTM data
+    """
+    path = pkg_resources.resource_filename("map2loop", "/_datasets/geodata_files/hamersley/dtm_rp.tif")
+    return gdal.Open(path)

map2loop/deformation_history.py

@@ -30,7 +30,8 @@ def __init__(self):
         """
         The initialiser for the deformation history. All attributes are defaulted
         """
-        self.minimum_fault_length_to_export = 500.0
+
+        self.minimum_fault_length_to_export = 0.0
         self.history = []
         self.fault_fault_relationships = []
 

map2loop/map2model_wrapper.py

@@ -133,9 +133,9 @@ def run(self, verbose_level: VerboseLevel = None):
         }
         # TODO: Simplify. Note: this is external so have to match fix to map2model module
         run_log = map2model.run(
-            os.path.join(self.map_data.tmp_path, "map2model_data"),
-            os.path.join(self.map_data.tmp_path, "map2model_data", "geology_wkt.csv"),
-            os.path.join(self.map_data.tmp_path, "map2model_data", "faults_wkt.csv"),
+            os.path.join(self.map_data.map2model_tmp_path),
+            os.path.join(self.map_data.map2model_tmp_path, "geology_wkt.csv"),
+            os.path.join(self.map_data.map2model_tmp_path, "faults_wkt.csv"),
             "",
             self.map_data.get_bounding_box(),
             map2model_code_map,
@@ -150,7 +150,7 @@ def run(self, verbose_level: VerboseLevel = None):
 
         # Parse units sorted
         units_sorted = pandas.read_csv(
-            os.path.join(self.map_data.tmp_path, "map2model_data", "units_sorted.txt"),
+            os.path.join(self.map_data.map2model_tmp_path, "units_sorted.txt"),
             header=None,
             sep=' ',
         )
@@ -162,7 +162,7 @@ def run(self, verbose_level: VerboseLevel = None):
         # Parse fault intersections
         out = []
         fault_fault_intersection_filename = os.path.join(
-            self.map_data.tmp_path, "map2model_data", "fault-fault-intersection.txt"
+            self.map_data.map2model_tmp_path, "fault-fault-intersection.txt"
         )
         if (
             os.path.isfile(fault_fault_intersection_filename)
@@ -189,7 +189,7 @@ def run(self, verbose_level: VerboseLevel = None):
         # Parse unit fault relationships
         out = []
         unit_fault_intersection_filename = os.path.join(
-            self.map_data.tmp_path, "map2model_data", "unit-fault-intersection.txt"
+            self.map_data.map2model_tmp_path, "unit-fault-intersection.txt"
         )
         if (
             os.path.isfile(unit_fault_intersection_filename)
@@ -212,11 +212,11 @@ def run(self, verbose_level: VerboseLevel = None):
         units = []
         links = []
         graph_filename = os.path.join(
-            self.map_data.tmp_path, "map2model_data", "graph_all_None.gml.txt"
+            self.map_data.map2model_tmp_path, "graph_all_None.gml.txt"
         )
         if os.path.isfile(graph_filename) and os.path.getsize(graph_filename) > 0:
             with open(
-                os.path.join(self.map_data.tmp_path, "map2model_data", "graph_all_None.gml.txt")
+                os.path.join(self.map_data.map2model_tmp_path, "graph_all_None.gml.txt")
             ) as file:
                 contents = file.read()
                 segments = contents.split("\n\n")

map2loop/mapdata.py

@@ -19,6 +19,7 @@
 import os
 from io import BytesIO
 from typing import Union
+import tempfile
 
 class MapData:
     """
@@ -52,15 +53,13 @@ class MapData:
         The filename of the json config file
     colour_filename: str
         The filename of the csv colour table file (columns are unit name and colour in #000000 form)
-    tmp_path: str
-        The path to the directory holding the temporary files
     verbose_level: m2l_enums.VerboseLevel
         A selection that defines how much console logging is output
     config: Config
         A link to the config structure which is defined in config.py
     """
 
-    def __init__(self, tmp_path: str = "", verbose_level: VerboseLevel = VerboseLevel.ALL):
+    def __init__(self, verbose_level: VerboseLevel = VerboseLevel.ALL):
         """
         The initialiser for the map data
 
@@ -76,7 +75,6 @@ def __init__(self, tmp_path: str = "", verbose_level: VerboseLevel = VerboseLeve
         self.basal_contacts = None
         self.sampled_contacts = None
         self.filenames = [None] * len(Datatype)
-        # self.output_filenames = [None] * len(Datatype)
         self.dirtyflags = [True] * len(Datatype)
         self.data_states = [Datastate.UNNAMED] * len(Datatype)
         self.working_projection = None
@@ -85,7 +83,6 @@ def __init__(self, tmp_path: str = "", verbose_level: VerboseLevel = VerboseLeve
         self.bounding_box_str = None
         self.config_filename = None
         self.colour_filename = None
-        self.tmp_path = tmp_path
         self.verbose_level = verbose_level
 
         self.config = Config()
@@ -99,7 +96,7 @@ def set_working_projection(self, projection):
                 The projection to use for map reprojection
         """
         if issubclass(type(projection), int):
-            projection = "EPSG:" + str(projection)
+            projection = f"EPSG:{str(projection)}"
             self.working_projection = projection
         elif issubclass(type(projection), str):
             self.working_projection = projection
@@ -473,17 +470,8 @@ def open_http_query(url: str):
             else:
                 return response
         except urllib.URLError:
-            return None
-
-    @beartype.beartype
-    def __check_and_create_tmp_path(self):
-        """
-        Create the temporary files directory if it is not valid
-        """
-        if not os.path.isdir(self.tmp_path):
-            os.mkdir(self.tmp_path)
-
-
+            return None        
+
 
     @beartype.beartype
     def __retrieve_tif(self, filename: str):
@@ -497,38 +485,38 @@ def __retrieve_tif(self, filename: str):
         Returns:
             _type_: The open geotiff in a gdal handler
         """
-        self.__check_and_create_tmp_path()
-
+
         # For gdal debugging use exceptions
         gdal.UseExceptions()
         bb_ll = tuple(float(coord) for coord in self.bounding_box_polygon.to_crs("EPSG:4326").geometry.total_bounds)
 
-        if filename.lower() == "aus" or filename.lower() == "au":
+        if filename.lower() in {"aus", "au"}:
 
             url = "http://services.ga.gov.au/gis/services/DEM_SRTM_1Second_over_Bathymetry_Topography/MapServer/WCSServer?"
             wcs = WebCoverageService(url, version="1.0.0")
 
             coverage = wcs.getCoverage(
                 identifier="1", bbox=bb_ll, format="GeoTIFF", crs=4326, width=2048, height=2048
             )
-            
+
             # This is stupid that gdal cannot read a byte stream and has to have a
             # file on the local system to open or otherwise create a gdal file
             # from scratch with Create
-
-            tmp_file = os.path.join(self.tmp_path, "StupidGDALLocalFile.tif")
-
+            import pathlib
+            tmp_file = pathlib.Path(tempfile.mkdtemp()) / pathlib.Path("temp.tif")
+            
             with open(tmp_file, "wb") as fh:
                 fh.write(coverage.read())
+
             tif = gdal.Open(tmp_file)
-        
+
         elif filename == "hawaii":
             import netCDF4
 
             bbox_str = (
                 f"[({str(bb_ll[1])}):1:({str(bb_ll[3])})][({str(bb_ll[0])}):1:({str(bb_ll[2])})]"
             )
-            
+
             filename = f"https://pae-paha.pacioos.hawaii.edu/erddap/griddap/srtm30plus_v11_land.nc?elev{bbox_str}"
             f = urllib.request.urlopen(filename)
             ds = netCDF4.Dataset("in-mem-file", mode="r", memory=f.read())
@@ -909,6 +897,7 @@ def parse_fault_map(self) -> tuple:
             tuple: A tuple of (bool: success/fail, str: failure message)
         """
         # Check type of loaded fault map
+
         if (
             self.raw_data[Datatype.FAULT] is None
             or type(self.raw_data[Datatype.FAULT]) is not geopandas.GeoDataFrame
@@ -917,6 +906,20 @@ def parse_fault_map(self) -> tuple:
 
         # Create new geodataframe
         faults = geopandas.GeoDataFrame(self.raw_data[Datatype.FAULT]["geometry"])
+
+        # crop by minimum fault length
+        self.minimum_fault_length = 0.0
+
+        try:
+            self.minimum_fault_length = self.config.fault_config["minimum_fault_length"]
+        except KeyError: ### if argument missing from config dictionary
+            print("Minimum fault length not defined in config dictionary; including all faults in map2loop project\n") ## TODO: pass this warning to logger
+
+
+        # Filter faults based on the minimum fault length
+        faults = faults.loc[faults.geometry.length > self.minimum_fault_length]
+
+
         config = self.config.fault_config
 
         if config["structtype_column"] in self.raw_data[Datatype.FAULT]:
@@ -1012,6 +1015,7 @@ def parse_fault_map(self) -> tuple:
         self.data[Datatype.FAULT] = faults
         return (False, "")
 
+
     @beartype.beartype
     def parse_fold_map(self) -> tuple:
         """
@@ -1121,15 +1125,15 @@ def save_raw_map_data(self, output_dir: str, datatype: Datatype, extension: str
                 The extension to use for the data. Defaults to ".csv".
         """
         try:
-            filename = os.path.join(output_dir, datatype.name + extension)
+            filename = pathlib.Path(output_dir) / f"{datatype.name}{extension}"
+
             if extension == ".csv":
-                # TODO: Add geopandas to pandas converter and then write csv
-                # self.raw_data[datatype].write_csv(filename)
-                print("GeoDataFrame to CSV conversion not implimented")
+                df = self.raw_data[datatype]
+                df.to_csv(filename, index=False)
             else:
                 self.raw_data[datatype].to_file(filename)
-        except Exception:
-            print(f"Failed to save {datatype.name} to file named {filename}\n")
+        except Exception as e:
+            print(f"Failed to save {datatype.name} to file named {filename}\nError: {str(e)}")
 
     @beartype.beartype
     def get_raw_map_data(self, datatype: Datatype):
@@ -1241,10 +1245,8 @@ def export_wkt_format_files(self):
         """
         # TODO: - Move away from tab seperators entirely (topology and map2model)
 
-        self.__check_and_create_tmp_path()
-        if not os.path.isdir(os.path.join(self.tmp_path, "map2model_data")):
-            os.mkdir(os.path.join(self.tmp_path, "map2model_data"))
-
+        self.map2model_tmp_path = pathlib.Path(tempfile.mkdtemp()) 
+
         # Check geology data status and export to a WKT format file
         self.load_map_data(Datatype.GEOLOGY)
         if type(self.data[Datatype.GEOLOGY]) is not geopandas.GeoDataFrame:
@@ -1277,7 +1279,7 @@ def export_wkt_format_files(self):
             geology["ROCKTYPE1"] = geology["ROCKTYPE1"].replace("", "None")
             geology["ROCKTYPE2"] = geology["ROCKTYPE2"].replace("", "None")
             geology.to_csv(
-                os.path.join(self.tmp_path, "map2model_data", "geology_wkt.csv"),
+                os.path.join(self.map2model_tmp_path, "geology_wkt.csv"),
                 sep="\t",
                 index=False,
             )
@@ -1294,7 +1296,7 @@ def export_wkt_format_files(self):
             faults = self.get_map_data(Datatype.FAULT).copy()
             faults.rename(columns={"geometry": "WKT"}, inplace=True)
             faults.to_csv(
-                os.path.join(self.tmp_path, "map2model_data", "faults_wkt.csv"),
+                os.path.join(self.map2model_tmp_path, "faults_wkt.csv"),
                 sep="\t",
                 index=False,
             )
@@ -1428,7 +1430,7 @@ def extract_all_contacts(self, save_contacts=True):
         if save_contacts:
             self.contacts = contacts
         return contacts
-
+    
     @beartype.beartype
     def extract_basal_contacts(self, stratigraphic_column: list, save_contacts=True):
         """
@@ -1440,19 +1442,46 @@ def extract_basal_contacts(self, stratigraphic_column: list, save_contacts=True)
         """
         units = stratigraphic_column
         basal_contacts = self.contacts.copy()
+
+        # check if the units in the strati colum are in the geology dataset, so that basal contacts can be built
+        # if not, stop the project
+        if any(
+            unit not in units for unit in basal_contacts["UNITNAME_1"].unique()
+        ):
+            missing_units = basal_contacts[~basal_contacts["UNITNAME_1"].isin(units)]["UNITNAME_1"].unique().tolist()
+            raise ValueError(
+                "There are units in stratigraphic column, but not in the Geology dataset: " +
+                ", ".join(missing_units) +
+                ". Please readjust the stratigraphic column if this is a user defined column."
+            )
+
+        # apply minimum lithological id between the two units 
         basal_contacts["ID"] = basal_contacts.apply(
             lambda row: min(units.index(row["UNITNAME_1"]), units.index(row["UNITNAME_2"])), axis=1
         )
+        # match the name of the unit with the minimum id
         basal_contacts["basal_unit"] = basal_contacts.apply(lambda row: units[row["ID"]], axis=1)
-        basal_contacts["distance"] = basal_contacts.apply(
+        # how many units apart are the two units? 
+        basal_contacts["stratigraphic_distance"] = basal_contacts.apply(
             lambda row: abs(units.index(row["UNITNAME_1"]) - units.index(row["UNITNAME_2"])), axis=1
         )
+        # if the units are more than 1 unit apart, the contact is abnormal (meaning that there is one (or more) unit(s) missing in between the two)
         basal_contacts["type"] = basal_contacts.apply(
-            lambda row: "ABNORMAL" if abs(row["distance"]) > 1 else "BASAL", axis=1
+            lambda row: "ABNORMAL" if abs(row["stratigraphic_distance"]) > 1 else "BASAL", axis=1
         )
+
         basal_contacts = basal_contacts[["ID", "basal_unit", "type", "geometry"]]
+
+        # added code to make sure that multi-line that touch each other are snapped and merged. 
+        # necessary for the reconstruction based on featureId
+        basal_contacts["geometry"] = [shapely.line_merge(shapely.snap(geo, geo, 1)) for geo in basal_contacts["geometry"]]
+
         if save_contacts:
-            self.basal_contacts = basal_contacts
+            #keep abnormal contacts as all_basal_contacts
+            self.all_basal_contacts = basal_contacts
+            #remove the abnormal contacts from basal contacts
+            self.basal_contacts = basal_contacts[basal_contacts["type"] == "BASAL"]
+
         return basal_contacts
 
     @beartype.beartype
@@ -1484,7 +1513,7 @@ def colour_units(
                 self.colour_filename = None
 
         if self.colour_filename is None:
-            print("No colour configuration file found. Assigning random colors to units")
+            print("\nNo colour configuration file found. Assigning random colors to units")
             missing_colour_n = len(stratigraphic_units["colour"])
             stratigraphic_units.loc[stratigraphic_units["colour"].isna(), "colour"] = (
                 generate_random_hex_colors(missing_colour_n)