sampler

Author: Noelle Cheng

map2loop/mapdata.py

@@ -1433,63 +1433,6 @@ def get_value_from_raster(self, datatype: Datatype, x, y):
         val = data.ReadAsArray(px, py, 1, 1)[0][0]
         return val
 
-    @beartype.beartype
-    def __value_from_raster(self, inv_geotransform, data, x: float, y: float):
-        """
-        Get the value from a raster dataset at the specified point
-
-        Args:
-            inv_geotransform (gdal.GeoTransform):
-                The inverse of the data's geotransform
-            data (numpy.array):
-                The raster data
-            x (float):
-                The easting coordinate of the value
-            y (float):
-                The northing coordinate of the value
-
-        Returns:
-            float or int: The value at the point specified
-        """
-        px = int(inv_geotransform[0] + inv_geotransform[1] * x + inv_geotransform[2] * y)
-        py = int(inv_geotransform[3] + inv_geotransform[4] * x + inv_geotransform[5] * y)
-        # Clamp values to the edges of raster if past boundary, similiar to GL_CLIP
-        px = max(px, 0)
-        px = min(px, data.shape[0] - 1)
-        py = max(py, 0)
-        py = min(py, data.shape[1] - 1)
-        return data[px][py]
-
-    @beartype.beartype
-    def get_value_from_raster_df(self, datatype: Datatype, df: pandas.DataFrame):
-        """
-        Add a 'Z' column to a dataframe with the heights from the 'X' and 'Y' coordinates
-
-        Args:
-            datatype (Datatype):
-                The datatype of the raster map to retrieve from
-            df (pandas.DataFrame):
-                The original dataframe with 'X' and 'Y' columns
-
-        Returns:
-            pandas.DataFrame: The modified dataframe
-        """
-        if len(df) <= 0:
-            df["Z"] = []
-            return df
-        data = self.get_map_data(datatype)
-        if data is None:
-            logger.warning("Cannot get value from data as data is not loaded")
-            return None
-
-        inv_geotransform = gdal.InvGeoTransform(data.GetGeoTransform())
-        data_array = numpy.array(data.GetRasterBand(1).ReadAsArray().T)
-
-        df["Z"] = df.apply(
-            lambda row: self.__value_from_raster(inv_geotransform, data_array, row["X"], row["Y"]),
-            axis=1,
-        )
-        return df
 
     @beartype.beartype
     def extract_all_contacts(self, save_contacts=True):

map2loop/project.py

@@ -1,6 +1,6 @@
 # internal imports
 from map2loop.fault_orientation import FaultOrientationNearest
-from .utils import hex_to_rgb
+from .utils import hex_to_rgb, set_z_values_from_raster_df
 from .m2l_enums import VerboseLevel, ErrorState, Datatype
 from .mapdata import MapData
 from .sampler import Sampler, SamplerDecimator, SamplerSpacing
@@ -506,23 +506,19 @@ def sample_map_data(self):
         logger.info(
             f"Sampling geology map data using {self.samplers[Datatype.GEOLOGY].sampler_label}"
         )
-        self.geology_samples = self.samplers[Datatype.GEOLOGY].sample(
-            self.map_data.get_map_data(Datatype.GEOLOGY), self.map_data
-        )
-        logger.info(
-            f"Sampling structure map data using {self.samplers[Datatype.STRUCTURE].sampler_label}"
-        )
-        self.structure_samples = self.samplers[Datatype.STRUCTURE].sample(
-            self.map_data.get_map_data(Datatype.STRUCTURE), self.map_data
-        )
+        geology_data = self.map_data.get_map_data(Datatype.GEOLOGY)
+        dtm_data = self.map_data.get_map_data(Datatype.DTM)
+
+        self.geology_samples = self.samplers[Datatype.GEOLOGY].sample(geology_data)
+        logger.info(f"Sampling structure map data using {self.samplers[Datatype.STRUCTURE].sampler_label}")
+
+        self.structure_samples = self.samplers[Datatype.STRUCTURE].sample(self.map_data.get_map_data(Datatype.STRUCTURE), dtm_data, geology_data)
         logger.info(f"Sampling fault map data using {self.samplers[Datatype.FAULT].sampler_label}")
-        self.fault_samples = self.samplers[Datatype.FAULT].sample(
-            self.map_data.get_map_data(Datatype.FAULT), self.map_data
-        )
+
+        self.fault_samples = self.samplers[Datatype.FAULT].sample(self.map_data.get_map_data(Datatype.FAULT))
         logger.info(f"Sampling fold map data using {self.samplers[Datatype.FOLD].sampler_label}")
-        self.fold_samples = self.samplers[Datatype.FOLD].sample(
-            self.map_data.get_map_data(Datatype.FOLD), self.map_data
-        )
+        
+        self.fold_samples = self.samplers[Datatype.FOLD].sample(self.map_data.get_map_data(Datatype.FOLD))
 
     def extract_geology_contacts(self):
         """
@@ -532,11 +528,9 @@ def extract_geology_contacts(self):
         self.map_data.extract_basal_contacts(self.stratigraphic_column.column)
 
         # sample the contacts
-        self.map_data.sampled_contacts = self.samplers[Datatype.GEOLOGY].sample(
-            self.map_data.basal_contacts
-        )
-
-        self.map_data.get_value_from_raster_df(Datatype.DTM, self.map_data.sampled_contacts)
+        self.map_data.sampled_contacts = self.samplers[Datatype.GEOLOGY].sample(self.map_data.basal_contacts)
+        dtm_data = self.map_data.get_map_data(Datatype.DTM)
+        set_z_values_from_raster_df(dtm_data, self.map_data.sampled_contacts)
 
     def calculate_stratigraphic_order(self, take_best=False):
         """
@@ -714,7 +708,8 @@ def calculate_fault_orientations(self):
                 self.map_data.get_map_data(Datatype.FAULT_ORIENTATION),
                 self.map_data,
             )
-            self.map_data.get_value_from_raster_df(Datatype.DTM, self.fault_orientations)
+            dtm_data = self.map_data.get_map_data(Datatype.DTM)
+            set_z_values_from_raster_df(dtm_data, self.fault_orientations)
         else:
             logger.warning(
                 "No fault orientation data found, skipping fault orientation calculation"
@@ -739,7 +734,8 @@ def summarise_fault_data(self):
         """
         Use the fault shapefile to make a summary of each fault by name
         """
-        self.map_data.get_value_from_raster_df(Datatype.DTM, self.fault_samples)
+        dtm_data = self.map_data.get_map_data(Datatype.DTM)
+        set_z_values_from_raster_df(dtm_data, self.fault_samples)
 
         self.deformation_history.summarise_data(self.fault_samples)
         self.deformation_history.faults = self.throw_calculator.compute(

map2loop/sampler.py

@@ -1,6 +1,7 @@
 # internal imports
 from .m2l_enums import Datatype
 from .mapdata import MapData
+from .utils import set_z_values_from_raster_df
 
 # external imports
 from abc import ABC, abstractmethod
@@ -10,6 +11,7 @@
 import shapely
 import numpy
 from typing import Optional
+from osgeo import gdal
 
 
 class Sampler(ABC):
@@ -38,7 +40,7 @@ def type(self):
     @beartype.beartype
     @abstractmethod
     def sample(
-        self, spatial_data: geopandas.GeoDataFrame, map_data: Optional[MapData] = None
+        self, spatial_data: geopandas.GeoDataFrame, dtm_data: Optional[geopandas.GeoDataFrame] = None, geology_data: Optional[geopandas.GeoDataFrame] = None
     ) -> pandas.DataFrame:
         """
         Execute sampling method (abstract method)
@@ -73,7 +75,7 @@ def __init__(self, decimation: int = 1):
 
     @beartype.beartype
     def sample(
-        self, spatial_data: geopandas.GeoDataFrame, map_data: Optional[MapData] = None
+        self, spatial_data: geopandas.GeoDataFrame, dtm_data: Optional[gdal.Dataset] = None, geology_data: Optional[geopandas.GeoDataFrame] = None
     ) -> pandas.DataFrame:
         """
         Execute sample method takes full point data, samples the data and returns the decimated points
@@ -87,10 +89,16 @@ def sample(
         data = spatial_data.copy()
         data["X"] = data.geometry.x
         data["Y"] = data.geometry.y
-        data["Z"] = map_data.get_value_from_raster_df(Datatype.DTM, data)["Z"]
-        data["layerID"] = geopandas.sjoin(
-            data, map_data.get_map_data(Datatype.GEOLOGY), how='left'
-        )['index_right']
+        if dtm_data is not None:
+            data["Z"] = set_z_values_from_raster_df(dtm_data, data)["Z"]
+        else:
+            data["Z"] = None
+        if geology_data is not None:
+            data["layerID"] = geopandas.sjoin(
+                data, geology_data, how='left'
+            )['index_right']
+        else:
+            data["layerID"] = None
         data.reset_index(drop=True, inplace=True)
 
         return pandas.DataFrame(data[:: self.decimation].drop(columns="geometry"))
@@ -118,7 +126,7 @@ def __init__(self, spacing: float = 50.0):
 
     @beartype.beartype
     def sample(
-        self, spatial_data: geopandas.GeoDataFrame, map_data: Optional[MapData] = None
+        self, spatial_data: geopandas.GeoDataFrame, dtm_data: Optional[geopandas.GeoDataFrame] = None, geology_data: Optional[geopandas.GeoDataFrame] = None
     ) -> pandas.DataFrame:
         """
         Execute sample method takes full point data, samples the data and returns the sampled points

map2loop/thickness_calculator.py

@@ -5,6 +5,7 @@
     calculate_endpoints,
     multiline_to_line,
     find_segment_strike_from_pt,
+    set_z_values_from_raster_df
 )
 from .m2l_enums import Datatype
 from .interpolators import DipDipDirectionInterpolator
@@ -271,7 +272,8 @@ def compute(
         # set the crs of the contacts to the crs of the units
         contacts = contacts.set_crs(crs=basal_contacts.crs)
         # get the elevation Z of the contacts
-        contacts = map_data.get_value_from_raster_df(Datatype.DTM, contacts)
+        dtm_data = map_data.get_map_data(Datatype.DTM)
+        contacts = set_z_values_from_raster_df(dtm_data, contacts)
         # update the geometry of the contact points to include the Z value
         contacts["geometry"] = contacts.apply(
             lambda row: shapely.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
@@ -299,7 +301,8 @@ def compute(
         # set the crs of the interpolated orientations to the crs of the units
         interpolated_orientations = interpolated_orientations.set_crs(crs=basal_contacts.crs)
         # get the elevation Z of the interpolated points
-        interpolated = map_data.get_value_from_raster_df(Datatype.DTM, interpolated_orientations)
+        dtm_data = map_data.get_map_data(Datatype.DTM)
+        interpolated = set_z_values_from_raster_df(dtm_data, interpolated_orientations)
         # update the geometry of the interpolated points to include the Z value
         interpolated["geometry"] = interpolated.apply(
             lambda row: shapely.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1

map2loop/utils.py

@@ -7,6 +7,7 @@
 import pandas
 import re
 import json
+from osgeo import gdal
 
 from .logging import getLogger
 logger = getLogger(__name__)
@@ -527,4 +528,62 @@ def update_from_legacy_file(
         with open(json_save_path, "w") as f:
             json.dump(parsed_data, f, indent=4)
 
-    return file_map
+    return file_map
+
+@beartype.beartype
+def value_from_raster(inv_geotransform, data, x: float, y: float):
+    """
+    Get the value from a raster dataset at the specified point
+
+    Args:
+        inv_geotransform (gdal.GeoTransform):
+            The inverse of the data's geotransform
+        data (numpy.array):
+            The raster data
+        x (float):
+            The easting coordinate of the value
+        y (float):
+            The northing coordinate of the value
+
+    Returns:
+        float or int: The value at the point specified
+    """
+    px = int(inv_geotransform[0] + inv_geotransform[1] * x + inv_geotransform[2] * y)
+    py = int(inv_geotransform[3] + inv_geotransform[4] * x + inv_geotransform[5] * y)
+    # Clamp values to the edges of raster if past boundary, similiar to GL_CLIP
+    px = max(px, 0)
+    px = min(px, data.shape[0] - 1)
+    py = max(py, 0)
+    py = min(py, data.shape[1] - 1)
+    return data[px][py]
+
+@beartype.beartype
+def set_z_values_from_raster_df(dtm_data: gdal.Dataset, df: pandas.DataFrame):
+    """
+    Add a 'Z' column to a dataframe with the heights from the 'X' and 'Y' coordinates
+
+    Args:
+        dtm_data (gdal.Dataset):
+            Dtm data from raster map
+        df (pandas.DataFrame):
+            The original dataframe with 'X' and 'Y' columns
+
+    Returns:
+        pandas.DataFrame: The modified dataframe
+    """
+    if len(df) <= 0:
+        df["Z"] = []
+        return df
+    
+    if dtm_data is None:
+        logger.warning("Cannot get value from data as data is not loaded")
+        return None
+
+    inv_geotransform = gdal.InvGeoTransform(dtm_data.GetGeoTransform())
+    data_array = numpy.array(dtm_data.GetRasterBand(1).ReadAsArray().T)
+
+    df["Z"] = df.apply(
+        lambda row: value_from_raster(inv_geotransform, data_array, row["X"], row["Y"]),
+        axis=1,
+    )
+    return df

pyproject.toml

@@ -6,7 +6,7 @@ name = 'map2loop'
 description = 'Generate 3D model data from 2D maps.'
 authors = [{name = 'Loop team'}]
 readme = 'README.md'
-requires-python = '>=3.8'
+requires-python = '>=3.8,<3.13'
 keywords = [ "earth sciences",
     "geology",
     "3-D modelling",