diff --git a/.github/workflows/conda.yml b/.github/workflows/conda.yml
index 7227a9b2..d60fc74e 100644
--- a/.github/workflows/conda.yml
+++ b/.github/workflows/conda.yml
@@ -33,7 +33,7 @@ jobs:
           conda install -c conda-forge conda-build scikit-build-core numpy anaconda-client conda-libmamba-solver -y
           conda config --set solver libmamba
           conda build --output-folder conda conda  --python ${{matrix.python-version}}
-          anaconda upload --label main conda/*/*.conda
+          anaconda upload --label main conda/*/*.tar.bz2
           
       - name: upload artifacts
         uses: actions/upload-artifact@v4
diff --git a/.github/workflows/linting_and_testing.yml b/.github/workflows/linting_and_testing.yml
index c432d6ce..9570a3bb 100644
--- a/.github/workflows/linting_and_testing.yml
+++ b/.github/workflows/linting_and_testing.yml
@@ -48,28 +48,23 @@ jobs:
         with:
             python-version: ${{ matrix.python-version }}
             conda-remove-defaults: "true"
-      
-      - name: Install dependencies for windows python 3.9 and 3.10
-        if: ${{ matrix.os == 'windows-latest' && (matrix.python-version == '3.9' || matrix.python-version == '3.10') }}
+
+
+      - name: Install dependencies for windows python 3.10
+        if: ${{ matrix.os == 'windows-latest' && matrix.python-version == '3.10' }}
         run: |
             conda run -n test conda info
             conda run -n test conda install -c loop3d -c conda-forge "gdal=3.4.3" python=${{ matrix.python-version }} -y
             conda run -n test conda install -c loop3d -c conda-forge --file dependencies.txt python=${{ matrix.python-version }} -y
             conda run -n test conda install pytest python=${{ matrix.python-version }} -y
 
-      - name: Install dependencies for windows python 3.11
-        if: ${{ matrix.os == 'windows-latest' && matrix.python-version == '3.11' }}
+      - name: Install dependencies for other environments
+        if: ${{ matrix.os != 'windows-latest' || matrix.python-version != '3.10' }}
         run: |
             conda run -n test conda info
-            conda run -n test conda install -c loop3d -c conda-forge gdal "netcdf4=*=nompi_py311*" python=${{ matrix.python-version }} -y
+            conda run -n test conda install -c loop3d -c conda-forge gdal python=${{ matrix.python-version }} -y
             conda run -n test conda install -c loop3d -c conda-forge --file dependencies.txt python=${{ matrix.python-version }} -y
             conda run -n test conda install pytest python=${{ matrix.python-version }} -y
-
-      - name: Install dependencies for other environments
-        if: ${{ !(matrix.os == 'windows-latest' && (matrix.python-version == '3.9' || matrix.python-version == '3.10' || matrix.python-version == '3.11')) }}
-        run: |
-            conda run -n test conda info
-            conda run -n test conda install -c loop3d -c conda-forge python=${{ matrix.python-version }} gdal pytest --file dependencies.txt -y
             
       - name: Install map2loop
         run: |
@@ -77,4 +72,4 @@ jobs:
 
       - name: Run tests
         run: |
-            conda run -n test  pytest -v
\ No newline at end of file
+            conda run -n test  pytest
\ No newline at end of file
diff --git a/map2loop/contact_extractor.py b/map2loop/contact_extractor.py
deleted file mode 100644
index 56bc104e..00000000
--- a/map2loop/contact_extractor.py
+++ /dev/null
@@ -1,96 +0,0 @@
-import geopandas
-import pandas
-import shapely
-from .logging import getLogger
-from typing import Optional
-
-logger = getLogger(__name__)
-
-class ContactExtractor:
-    def __init__(self, geology: geopandas.GeoDataFrame, faults: Optional[geopandas.GeoDataFrame] = None):
-        self.geology = geology
-        self.faults = faults
-        self.contacts = None
-        self.basal_contacts = None
-        self.all_basal_contacts = None
-
-    def extract_all_contacts(self, save_contacts: bool = True) -> geopandas.GeoDataFrame:
-        logger.info("Extracting contacts")
-        geology = self.geology.copy()
-        geology = geology.dissolve(by="UNITNAME", as_index=False)
-        geology = geology[~geology["INTRUSIVE"]]
-        geology = geology[~geology["SILL"]]
-        if self.faults is not None:
-            faults = self.faults.copy()
-            faults["geometry"] = faults.buffer(50)
-            geology = geopandas.overlay(geology, faults, how="difference", keep_geom_type=False)
-        units = geology["UNITNAME"].unique().tolist()
-        column_names = ["UNITNAME_1", "UNITNAME_2", "geometry"]
-        contacts = geopandas.GeoDataFrame(crs=geology.crs, columns=column_names, data=None)
-        while len(units) > 1:
-            unit1 = units[0]
-            units = units[1:]
-            for unit2 in units:
-                if unit1 != unit2:
-                    join = geopandas.overlay(
-                        geology[geology["UNITNAME"] == unit1],
-                        geology[geology["UNITNAME"] == unit2],
-                        keep_geom_type=False,
-                    )[column_names]
-                    join["geometry"] = join.buffer(1)
-                    buffered = geology[geology["UNITNAME"] == unit2][["geometry"]].copy()
-                    buffered["geometry"] = buffered.boundary
-                    end = geopandas.overlay(buffered, join, keep_geom_type=False)
-                    if len(end):
-                        contacts = pandas.concat([contacts, end], ignore_index=True)
-        contacts["length"] = [row.length for row in contacts["geometry"]]
-        if save_contacts:
-            self.contacts = contacts
-        return contacts
-
-    def extract_basal_contacts(self, 
-                               stratigraphic_column: list, 
-                               save_contacts: bool = True) -> geopandas.GeoDataFrame:
-        
-        logger.info("Extracting basal contacts")
-        units = stratigraphic_column
-        
-        if self.contacts is None:
-            self.extract_all_contacts(save_contacts=True)
-            basal_contacts = self.contacts.copy()
-        else:
-            basal_contacts = self.contacts.copy()
-        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()
-            )
-            logger.error(
-                "There are units in the Geology dataset, but not in the stratigraphic column: "
-                + ", ".join(missing_units)
-                + ". Please readjust the stratigraphic column if this is a user defined column."
-            )
-            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."
-            )
-        basal_contacts["ID"] = basal_contacts.apply(
-            lambda row: min(units.index(row["UNITNAME_1"]), units.index(row["UNITNAME_2"])), axis=1
-        )
-        basal_contacts["basal_unit"] = basal_contacts.apply(lambda row: units[row["ID"]], axis=1)
-        basal_contacts["stratigraphic_distance"] = basal_contacts.apply(
-            lambda row: abs(units.index(row["UNITNAME_1"]) - units.index(row["UNITNAME_2"])), axis=1
-        )
-        basal_contacts["type"] = basal_contacts.apply(
-            lambda row: "ABNORMAL" if abs(row["stratigraphic_distance"]) > 1 else "BASAL", axis=1
-        )
-        basal_contacts = basal_contacts[["ID", "basal_unit", "type", "geometry"]]
-        basal_contacts["geometry"] = [
-            shapely.line_merge(shapely.snap(geo, geo, 1)) for geo in basal_contacts["geometry"]
-        ]
-        if save_contacts:
-            self.all_basal_contacts = basal_contacts
-            self.basal_contacts = basal_contacts[basal_contacts["type"] == "BASAL"]
-        return basal_contacts
diff --git a/map2loop/data_checks.py b/map2loop/data_checks.py
index 5d766267..732464d8 100644
--- a/map2loop/data_checks.py
+++ b/map2loop/data_checks.py
@@ -471,13 +471,11 @@ def check_keys(d: dict, parent_key=""):
         "fault": {
             "structtype_column",  "fault_text",  "dip_null_value",
             "dipdir_flag", "dipdir_column",  "dip_column",  "orientation_type",
-            "dipestimate_column",  "dipestimate_text","displacement_column", 
-            "displacement_text",  "name_column","objectid_column", "featureid_column", "minimum_fault_length",
-            "fault_length_column", "fault_length_text", "ignore_fault_codes",
+            "dipestimate_column",  "dipestimate_text",  "name_column",
+            "objectid_column", "minimum_fault_length", "ignore_fault_codes",
         },
         "fold": {
-            "structtype_column", "fold_text", "axial_plane_dipdir_column", 
-            "axial_plane_dip_column","tightness_column", "description_column",
+            "structtype_column", "fold_text", "description_column",
             "synform_text", "foldname_column","objectid_column",
         },
     }
@@ -855,7 +853,6 @@ def validate_structtype_column(
     if text_keys:
         for text_key, config_key in text_keys.items():
             text_value = config.get(config_key, None)
-            text_pattern = '|'.join(text_value) if text_value else None
             if text_value:
                 if not isinstance(text_value, str):
                     error_msg = (
@@ -865,7 +862,7 @@ def validate_structtype_column(
                     logger.error(error_msg)
                     return (True, error_msg)
                 
-                if not geodata[structtype_column].str.contains(text_pattern, case=False, regex=True, na=False).any():
+                if not geodata[structtype_column].str.contains(text_value, na=False).any():
                     if text_key == "synform_text":
                         warning_msg = (
                             f"Datatype {datatype_name.upper()}: The '{text_key}' value '{text_value}' is not found in column '{structtype_column}'. "
diff --git a/map2loop/topology.py b/map2loop/map2model_wrapper.py
similarity index 72%
rename from map2loop/topology.py
rename to map2loop/map2model_wrapper.py
index d7a90dc4..115b8702 100644
--- a/map2loop/topology.py
+++ b/map2loop/map2model_wrapper.py
@@ -1,26 +1,24 @@
-from typing import Optional
 # internal imports
 from .m2l_enums import VerboseLevel
-from .contact_extractor import ContactExtractor
 
 # external imports
-import geopandas
-import pandas
-import numpy
+import geopandas as gpd
+import pandas as pd
+import numpy as np
 
 from .logging import getLogger
 
 logger = getLogger(__name__)
 
 
-class Topology:
+class Map2ModelWrapper:
     """
-    Topology calculator for units and faults relationships
+    A wrapper around map2model functionality
 
     Attributes
     ----------
     sorted_units: None or list
-        the stratigraphic column
+        map2model's estimate of the stratigraphic column
     fault_fault_relationships: None or pandas.DataFrame
         data frame of fault to fault relationships with columns ["Fault1", "Fault2", "Type", "Angle"]
     unit_fault_relationships: None or pandas.DataFrame
@@ -34,9 +32,7 @@ class Topology:
     """
 
     def __init__(
-        self, geology_data: geopandas.GeoDataFrame, 
-        fault_data: Optional[geopandas.GeoDataFrame] = None,  
-        verbose_level: VerboseLevel = VerboseLevel.NONE
+        self, map_data, *,  verbose_level: VerboseLevel = VerboseLevel.NONE
     ):
         """
         The initialiser for the map2model wrapper
@@ -51,8 +47,7 @@ def __init__(
         self._fault_fault_relationships = None
         self._unit_fault_relationships = None
         self._unit_unit_relationships = None
-        self.geology_data = geology_data
-        self.fault_data = fault_data
+        self.map_data = map_data
         self.verbose_level = verbose_level
         self.buffer_radius = 500
 
@@ -79,13 +74,13 @@ def unit_unit_relationships(self):
 
     def reset(self):
         """
-        Reset before the next run
+        Reset the wrapper to before the map2model process
         """
-        logger.info("Resetting topology calculator")
+        logger.info("Resetting map2model wrapper")
         self.sorted_units = None
-        self._fault_fault_relationships = None
-        self._unit_fault_relationships = None
-        self._unit_unit_relationships = None
+        self.fault_fault_relationships = None
+        self.unit_fault_relationships = None
+        self.unit_unit_relationships = None
 
     def get_sorted_units(self):
         """
@@ -129,23 +124,23 @@ def get_unit_unit_relationships(self):
 
     def _calculate_fault_fault_relationships(self):
 
-        faults = self.fault_data.copy()
+        faults = self.map_data.FAULT.copy()
         # reset index so that we can index the adjacency matrix with the index
         faults.reset_index(inplace=True)
         buffers = faults.buffer(self.buffer_radius)
         # create the adjacency matrix
-        intersection = geopandas.sjoin(
-            geopandas.GeoDataFrame(geometry=buffers), geopandas.GeoDataFrame(geometry=faults["geometry"])
+        intersection = gpd.sjoin(
+            gpd.GeoDataFrame(geometry=buffers), gpd.GeoDataFrame(geometry=faults["geometry"])
         )
         intersection["index_left"] = intersection.index
         intersection.reset_index(inplace=True)
 
-        adjacency_matrix = numpy.zeros((faults.shape[0], faults.shape[0]), dtype=bool)
+        adjacency_matrix = np.zeros((faults.shape[0], faults.shape[0]), dtype=bool)
         adjacency_matrix[
             intersection.loc[:, "index_left"], intersection.loc[:, "index_right"]
         ] = True
-        f1, f2 = numpy.where(numpy.tril(adjacency_matrix, k=-1))
-        df = pandas.DataFrame(
+        f1, f2 = np.where(np.tril(adjacency_matrix, k=-1))
+        df = pd.DataFrame(
             {'Fault1': faults.loc[f1, 'ID'].to_list(), 'Fault2': faults.loc[f2, 'ID'].to_list()}
         )
         df['Angle'] = 60  # make it big to prevent LS from making splays
@@ -156,29 +151,26 @@ def _calculate_fault_unit_relationships(self):
         """Calculate unit/fault relationships using geopandas sjoin.
         This will return
         """
-        units = self.geology_data["UNITNAME"].unique()
-        faults = self.fault_data.copy().reset_index().drop(columns=['index'])
-        adjacency_matrix = numpy.zeros((len(units), faults.shape[0]), dtype=bool)
+        units = self.map_data.GEOLOGY["UNITNAME"].unique()
+        faults = self.map_data.FAULT.copy().reset_index().drop(columns=['index'])
+        adjacency_matrix = np.zeros((len(units), faults.shape[0]), dtype=bool)
         for i, u in enumerate(units):
-            unit = self.geology_data[self.geology_data["UNITNAME"] == u]
-            intersection = geopandas.sjoin(
-                geopandas.GeoDataFrame(geometry=faults["geometry"]),
-                geopandas.GeoDataFrame(geometry=unit["geometry"]),
+            unit = self.map_data.GEOLOGY[self.map_data.GEOLOGY["UNITNAME"] == u]
+            intersection = gpd.sjoin(
+                gpd.GeoDataFrame(geometry=faults["geometry"]),
+                gpd.GeoDataFrame(geometry=unit["geometry"]),
             )
             intersection["index_left"] = intersection.index
             intersection.reset_index(inplace=True)
             adjacency_matrix[i, intersection.loc[:, "index_left"]] = True
-        u, f = numpy.where(adjacency_matrix)
-        df = pandas.DataFrame({"Unit": units[u].tolist(), "Fault": faults.loc[f, "ID"].to_list()})
+        u, f = np.where(adjacency_matrix)
+        df = pd.DataFrame({"Unit": units[u].tolist(), "Fault": faults.loc[f, "ID"].to_list()})
         self._unit_fault_relationships = df
 
     def _calculate_unit_unit_relationships(self):
-        extractor = ContactExtractor(
-            self.geology_data,
-            self.fault_data,
-        )
-        contacts = extractor.extract_all_contacts()
-        self._unit_unit_relationships = contacts.copy().drop(
+        if self.map_data.contacts is None:
+            self.map_data.extract_all_contacts()
+        self._unit_unit_relationships = self.map_data.contacts.copy().drop(
             columns=['length', 'geometry']
         )
         return self._unit_unit_relationships
diff --git a/map2loop/mapdata.py b/map2loop/mapdata.py
index 5ae780c1..4137af27 100644
--- a/map2loop/mapdata.py
+++ b/map2loop/mapdata.py
@@ -596,11 +596,7 @@ def __retrieve_tif(self, filename: str):
             )
 
             filename = f"https://pae-paha.pacioos.hawaii.edu/erddap/griddap/srtm30plus_v11_land.nc?elev{bbox_str}"
-            try:
-                f = urllib.request.urlopen(filename)
-            except urllib.error.URLError:
-                logger.error(f"Failed to open remote file {filename}")
-                return None
+            f = urllib.request.urlopen(filename)
             ds = netCDF4.Dataset("in-mem-file", mode="r", memory=f.read())
             spatial = [
                 ds.geospatial_lon_min,
@@ -625,13 +621,7 @@ def __retrieve_tif(self, filename: str):
             tif.GetRasterBand(1).WriteArray(numpy.flipud(ds.variables["elev"][:][:]))
         elif filename.startswith("http"):
             logger.info(f'Opening remote file {filename}')
-            try:
-                image_data = self.open_http_query(filename)
-            except urllib.error.URLError:
-                logger.error(f"Failed to open remote file {filename}")
-                return None
-            if image_data is None:
-                return None
+            image_data = self.open_http_query(filename)
             mmap_name = f"/vsimem/{str(uuid4())}"
             gdal.FileFromMemBuffer(mmap_name, image_data.read())
             tif = gdal.Open(mmap_name)
@@ -655,9 +645,6 @@ def load_raster_map_data(self, datatype: Datatype):
             if self.data_states[datatype] == Datastate.UNLOADED:
                 # Load data from file
                 self.data[datatype] = self.__retrieve_tif(self.filenames[datatype])
-                if self.data[datatype] is None:
-                    logger.error(f"Failed to load raster data for {datatype.name}")
-                    return
                 self.data_states[datatype] = Datastate.LOADED
             if self.data_states[datatype] == Datastate.LOADED:
                 # Reproject raster to required CRS
@@ -672,7 +659,6 @@ def load_raster_map_data(self, datatype: Datatype):
                     )
                 except Exception:
                     logger.error(f"Warp failed for {datatype.name}\n")
-                    return
                 self.data_states[datatype] = Datastate.REPROJECTED
             if self.data_states[datatype] == Datastate.REPROJECTED:
                 # Clip raster image to bounding polygon
@@ -682,9 +668,6 @@ def load_raster_map_data(self, datatype: Datatype):
                     self.bounding_box["maxx"],
                     self.bounding_box["miny"],
                 ]
-                if self.data[datatype] is None:
-                    logger.error(f"No raster data available for {datatype.name}")
-                    return
                 self.data[datatype] = gdal.Translate(
                     "",
                     self.data[datatype],
@@ -1465,7 +1448,170 @@ 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):
+        """
+        Extract the contacts between units in the geology GeoDataFrame
+        """
+        logger.info("Extracting contacts")
+        geology = self.get_map_data(Datatype.GEOLOGY).copy()
+        geology = geology.dissolve(by="UNITNAME", as_index=False)
+        # Remove intrusions
+        geology = geology[~geology["INTRUSIVE"]]
+        geology = geology[~geology["SILL"]]
+        # Remove faults from contact geomety
+        if self.get_map_data(Datatype.FAULT) is not None:
+            faults = self.get_map_data(Datatype.FAULT).copy()
+            faults["geometry"] = faults.buffer(50)
+            geology = geopandas.overlay(geology, faults, how="difference", keep_geom_type=False)
+        units = geology["UNITNAME"].unique()
+        column_names = ["UNITNAME_1", "UNITNAME_2", "geometry"]
+        contacts = geopandas.GeoDataFrame(crs=geology.crs, columns=column_names, data=None)
+        while len(units) > 1:
+            unit1 = units[0]
+            units = units[1:]
+            for unit2 in units:
+                if unit1 != unit2:
+                    # print(f'contact: {unit1} and {unit2}')
+                    join = geopandas.overlay(
+                        geology[geology["UNITNAME"] == unit1],
+                        geology[geology["UNITNAME"] == unit2],
+                        keep_geom_type=False,
+                    )[column_names]
+                    join["geometry"] = join.buffer(1)
+                    buffered = geology[geology["UNITNAME"] == unit2][["geometry"]].copy()
+                    buffered["geometry"] = buffered.boundary
+                    end = geopandas.overlay(buffered, join, keep_geom_type=False)
+                    if len(end):
+                        contacts = pandas.concat([contacts, end], ignore_index=True)
+        # contacts["TYPE"] = "UNKNOWN"
+        contacts["length"] = [row.length for row in contacts["geometry"]]
+        # print('finished extracting contacts')
+        if save_contacts:
+            self.contacts = contacts
+        return contacts
+
+    @beartype.beartype
+    def extract_basal_contacts(self, stratigraphic_column: list, save_contacts=True):
+        """
+        Identify the basal unit of the contacts based on the stratigraphic column
+
+        Args:
+            stratigraphic_column (list):
+                The stratigraphic column to use
+        """
+        logger.info("Extracting basal contacts")
+        
+        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()
+            )
+            logger.error(
+                "There are units in the Geology dataset, but not in the stratigraphic column: "
+                + ", ".join(missing_units)
+                + ". Please readjust the stratigraphic column if this is a user defined column."
+            )
+            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)
+        # 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["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:
+            # 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
     def colour_units(
diff --git a/map2loop/project.py b/map2loop/project.py
index 2e1b46e6..d9cfbb83 100644
--- a/map2loop/project.py
+++ b/map2loop/project.py
@@ -1,9 +1,8 @@
 # internal imports
 from map2loop.fault_orientation import FaultOrientationNearest
-from .utils import hex_to_rgb, set_z_values_from_raster_df
+from .utils import hex_to_rgb
 from .m2l_enums import VerboseLevel, ErrorState, Datatype
 from .mapdata import MapData
-from .contact_extractor import ContactExtractor
 from .sampler import Sampler, SamplerDecimator, SamplerSpacing
 from .thickness_calculator import InterpolatedStructure, ThicknessCalculator
 from .throw_calculator import ThrowCalculator, ThrowCalculatorAlpha
@@ -11,7 +10,7 @@
 from .sorter import Sorter, SorterAgeBased, SorterAlpha, SorterUseNetworkX, SorterUseHint
 from .stratigraphic_column import StratigraphicColumn
 from .deformation_history import DeformationHistory
-from .topology import Topology
+from .map2model_wrapper import Map2ModelWrapper
 from .data_checks import validate_config_dictionary
 
 # external imports
@@ -41,14 +40,14 @@ class Project(object):
     verbose_level: m2l_enums.VerboseLevel
         A selection that defines how much console logging is output
     samplers: Sampler
-        A list of samplers used to extract point samples from polygons or line segments. Indexed by m2l_enum
+        A list of samplers used to extract point samples from polyonal or line segments. Indexed by m2l_enum.Dataype
     sorter: Sorter
         The sorting algorithm to use for calculating the stratigraphic column
     loop_filename: str
         The name of the loop project file used in this project
     map_data: MapData
         The structure that holds all map and dtm data
-    map2model: Topology
+    map2model: Map2ModelWrapper
         A wrapper around the map2model module that extracts unit and fault adjacency
     stratigraphic_column: StratigraphicColumn
         The structure that holds the unit information and ordering
@@ -139,18 +138,18 @@ def __init__(
         self.samplers = [SamplerDecimator()] * len(Datatype)
         self.set_default_samplers()
         self.bounding_box = bounding_box
-        self.contact_extractor = None
         self.sorter = SorterUseHint()
+        self.thickness_calculator = [InterpolatedStructure()]
         self.throw_calculator = ThrowCalculatorAlpha()
         self.fault_orientation = FaultOrientationNearest()
         self.map_data = MapData(verbose_level=verbose_level)
+        self.map2model = Map2ModelWrapper(self.map_data)
         self.stratigraphic_column = StratigraphicColumn()
         self.deformation_history = DeformationHistory(project=self)
         self.loop_filename = loop_project_filename
         self.overwrite_lpf = overwrite_loopprojectfile
         self.active_thickness = None
         
-        
         # initialise the dataframes to store data in
         self.fault_orientations = pandas.DataFrame(
             columns=["ID", "DIPDIR", "DIP", "X", "Y", "Z", "featureId"]
@@ -240,14 +239,6 @@ def __init__(
         # Populate the stratigraphic column and deformation history from map data
         self.stratigraphic_column.populate(self.map_data.get_map_data(Datatype.GEOLOGY))
         self.deformation_history.populate(self.map_data.get_map_data(Datatype.FAULT))
-        self.topology = Topology(
-            self.map_data.get_map_data(Datatype.GEOLOGY), 
-            self.map_data.get_map_data(Datatype.FAULT)
-            )
-        self.thickness_calculator = [InterpolatedStructure(
-            dtm_data=self.map_data.get_map_data(Datatype.DTM),
-            bounding_box=self.bounding_box,
-        )]
 
 
     @beartype.beartype
@@ -512,52 +503,45 @@ def sample_map_data(self):
         """
         Use the samplers to extract points along polylines or unit boundaries
         """
-        geology_data = self.map_data.get_map_data(Datatype.GEOLOGY)
-        dtm_data = self.map_data.get_map_data(Datatype.DTM)
-        
-        logger.info(f"Sampling geology map data using {self.samplers[Datatype.GEOLOGY].sampler_label}")
-        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.samplers[Datatype.STRUCTURE].dtm_data = dtm_data
-        self.samplers[Datatype.STRUCTURE].geology_data = geology_data
-        self.structure_samples = self.samplers[Datatype.STRUCTURE].sample(self.map_data.get_map_data(Datatype.STRUCTURE))
-
+        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
+        )
         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.fault_samples = self.samplers[Datatype.FAULT].sample(
+            self.map_data.get_map_data(Datatype.FAULT), self.map_data
+        )
         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.fold_samples = self.samplers[Datatype.FOLD].sample(
+            self.map_data.get_map_data(Datatype.FOLD), self.map_data
+        )
 
     def extract_geology_contacts(self):
         """
         Use the stratigraphic column, and fault and geology data to extract points along contacts
         """
         # Use stratigraphic column to determine basal contacts
-        if self.contact_extractor is None:
-            self.contact_extractor = ContactExtractor(
-                self.map_data.get_map_data(Datatype.GEOLOGY),
-                self.map_data.get_map_data(Datatype.FAULT),
-            )
-            self.contact_extractor.extract_all_contacts()
-
-        basal_contacts = self.contact_extractor.extract_basal_contacts(self.stratigraphic_column.column)
+        self.map_data.extract_basal_contacts(self.stratigraphic_column.column)
 
         # sample the contacts
-        self.map_data.sampled_contacts = self.samplers[Datatype.GEOLOGY].sample(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)
+        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)
 
     def calculate_stratigraphic_order(self, take_best=False):
         """
         Use unit relationships, unit ages and the sorter to create a stratigraphic column
         """
-        if self.contact_extractor is None:
-            self.contact_extractor = ContactExtractor(
-                self.map_data.get_map_data(Datatype.GEOLOGY),
-                self.map_data.get_map_data(Datatype.FAULT),
-            )
-            self.contact_extractor.extract_all_contacts()
         if take_best:
             sorters = [SorterUseHint(), SorterAgeBased(), SorterAlpha(), SorterUseNetworkX()]
             logger.info(
@@ -567,16 +551,14 @@ def calculate_stratigraphic_order(self, take_best=False):
             columns = [
                 sorter.sort(
                     self.stratigraphic_column.stratigraphicUnits,
-                    self.topology.get_unit_unit_relationships(),
-                    self.contact_extractor.contacts,
+                    self.map2model.get_unit_unit_relationships(),
+                    self.map_data.contacts,
                     self.map_data,
                 )
                 for sorter in sorters
             ]
             basal_contacts = [
-                self.contact_extractor.extract_basal_contacts(
-                    column, save_contacts=False
-                )
+                self.map_data.extract_basal_contacts(column, save_contacts=False)
                 for column in columns
             ]
             basal_lengths = [
@@ -599,8 +581,8 @@ def calculate_stratigraphic_order(self, take_best=False):
             logger.info(f'Calculating stratigraphic column using sorter {self.sorter.sorter_label}')
             self.stratigraphic_column.column = self.sorter.sort(
                 self.stratigraphic_column.stratigraphicUnits,
-                self.topology.get_unit_unit_relationships(),
-                self.contact_extractor.contacts,
+                self.map2model.get_unit_unit_relationships(),
+                self.map_data.contacts,
                 self.map_data,
             )
 
@@ -667,7 +649,7 @@ def get_thickness_calculator(self) -> List[str]:
                 "self.thickness_calculator must be either a list of objects or a single object with a thickness_calculator_label attribute"
             )
 
-    def calculate_unit_thicknesses(self, basal_contacts):
+    def calculate_unit_thicknesses(self):
         """
         Calculates the unit thickness statistics (mean, median, standard deviation) for each stratigraphic unit
         in the stratigraphic column using the provided thickness calculators.
@@ -694,15 +676,13 @@ def calculate_unit_thicknesses(self, basal_contacts):
         labels = []
 
         for calculator in self.thickness_calculator:
-            calculator.dtm_data = self.map_data.get_map_data(Datatype.DTM)
-            calculator.bounding_box = self.bounding_box
+
             result = calculator.compute(
                 self.stratigraphic_column.stratigraphicUnits,
                 self.stratigraphic_column.column,
-                basal_contacts,
+                self.map_data.basal_contacts,
                 self.structure_samples,
-                self.map_data.get_map_data(Datatype.GEOLOGY),
-                self.map_data.sampled_contacts,
+                self.map_data,
             )[['ThicknessMean', 'ThicknessMedian', 'ThicknessStdDev']].to_numpy()
 
             label = calculator.thickness_calculator_label
@@ -734,8 +714,7 @@ def calculate_fault_orientations(self):
                 self.map_data.get_map_data(Datatype.FAULT_ORIENTATION),
                 self.map_data,
             )
-            dtm_data = self.map_data.get_map_data(Datatype.DTM)
-            set_z_values_from_raster_df(dtm_data, self.fault_orientations)
+            self.map_data.get_value_from_raster_df(Datatype.DTM, self.fault_orientations)
         else:
             logger.warning(
                 "No fault orientation data found, skipping fault orientation calculation"
@@ -760,14 +739,13 @@ def summarise_fault_data(self):
         """
         Use the fault shapefile to make a summary of each fault by name
         """
-        dtm_data = self.map_data.get_map_data(Datatype.DTM)
-        set_z_values_from_raster_df(dtm_data, self.fault_samples)
+        self.map_data.get_value_from_raster_df(Datatype.DTM, self.fault_samples)
 
         self.deformation_history.summarise_data(self.fault_samples)
         self.deformation_history.faults = self.throw_calculator.compute(
             self.deformation_history.faults,
             self.stratigraphic_column.column,
-            self.contact_extractor.basal_contacts,
+            self.map_data.basal_contacts,
             self.map_data,
         )
         logger.info(f'There are {self.deformation_history.faults.shape[0]} faults in the dataset')
@@ -785,16 +763,12 @@ def run_all(self, user_defined_stratigraphic_column=None, take_best=False):
             logger.info(f'User defined stratigraphic column: {user_defined_stratigraphic_column}')
 
         # Calculate contacts before stratigraphic column
-        self.contact_extractor = ContactExtractor(
-            self.map_data.get_map_data(Datatype.GEOLOGY),
-            self.map_data.get_map_data(Datatype.FAULT),
-        )
-        self.map_data.contacts = self.contact_extractor.extract_all_contacts()
+        self.map_data.extract_all_contacts()
 
         # Calculate the stratigraphic column
         if issubclass(type(user_defined_stratigraphic_column), list):
             self.stratigraphic_column.column = user_defined_stratigraphic_column
-            self.topology.run()  # if we use a user defined stratigraphic column, we still need to calculate the results of map2model
+            self.map2model.run()  # if we use a user defined stratigraphic column, we still need to calculate the results of map2model
         else:
             if user_defined_stratigraphic_column is not None:
                 logger.warning(
@@ -806,7 +780,7 @@ def run_all(self, user_defined_stratigraphic_column=None, take_best=False):
         # Calculate basal contacts based on stratigraphic column
         self.extract_geology_contacts()
         self.sample_map_data()
-        self.calculate_unit_thicknesses(self.contact_extractor.basal_contacts)
+        self.calculate_unit_thicknesses()
         self.calculate_fault_orientations()
         self.summarise_fault_data()
         self.apply_colour_to_units()
@@ -1038,9 +1012,9 @@ def save_into_projectfile(self):
         observations["dipPolarity"] = self.structure_samples["OVERTURNED"]
         LPF.Set(self.loop_filename, "stratigraphicObservations", data=observations)
 
-        if self.topology.fault_fault_relationships is not None:
+        if self.map2model.fault_fault_relationships is not None:
             ff_relationships = self.deformation_history.get_fault_relationships_with_ids(
-                self.topology.fault_fault_relationships
+                self.map2model.fault_fault_relationships
             )
             relationships = numpy.zeros(len(ff_relationships), LPF.eventRelationshipType)
 
@@ -1079,7 +1053,7 @@ def draw_geology_map(self, points: pandas.DataFrame = None, overlay: str = ""):
             base = geol.plot(color=geol["colour_rgba"])
         if overlay != "":
             if overlay == "basal_contacts":
-                self.contact_extractor.basal_contacts[self.contact_extractor.basal_contacts["type"] == "BASAL"].plot(
+                self.map_data.basal_contacts[self.map_data.basal_contacts["type"] == "BASAL"].plot(
                     ax=base
                 )
 
diff --git a/map2loop/sampler.py b/map2loop/sampler.py
index b4c7835c..01600566 100644
--- a/map2loop/sampler.py
+++ b/map2loop/sampler.py
@@ -1,5 +1,6 @@
 # internal imports
-from .utils import set_z_values_from_raster_df
+from .m2l_enums import Datatype
+from .mapdata import MapData
 
 # external imports
 from abc import ABC, abstractmethod
@@ -9,7 +10,6 @@
 import shapely
 import numpy
 from typing import Optional
-from osgeo import gdal
 
 
 class Sampler(ABC):
@@ -20,13 +20,11 @@ class Sampler(ABC):
         ABC (ABC): Derived from Abstract Base Class
     """
 
-    def __init__(self, dtm_data=None, geology_data=None):
+    def __init__(self):
         """
         Initialiser of for Sampler
         """
         self.sampler_label = "SamplerBaseClass"
-        self.dtm_data = dtm_data
-        self.geology_data = geology_data
 
     def type(self):
         """
@@ -40,7 +38,7 @@ def type(self):
     @beartype.beartype
     @abstractmethod
     def sample(
-        self, spatial_data: geopandas.GeoDataFrame
+        self, spatial_data: geopandas.GeoDataFrame, map_data: Optional[MapData] = None
     ) -> pandas.DataFrame:
         """
         Execute sampling method (abstract method)
@@ -62,24 +60,20 @@ class SamplerDecimator(Sampler):
     """
 
     @beartype.beartype
-    def __init__(self, decimation: int = 1, dtm_data: Optional[gdal.Dataset] = None, geology_data: Optional[geopandas.GeoDataFrame] = None):
+    def __init__(self, decimation: int = 1):
         """
         Initialiser for decimator sampler
 
         Args:
             decimation (int, optional): stride of the points to sample. Defaults to 1.
-            dtm_data (Optional[gdal.Dataset], optional): digital terrain map data. Defaults to None.
-            geology_data (Optional[geopandas.GeoDataFrame], optional): geology data. Defaults to None.
         """
-        super().__init__(dtm_data, geology_data)
         self.sampler_label = "SamplerDecimator"
         decimation = max(decimation, 1)
         self.decimation = decimation
 
-
     @beartype.beartype
     def sample(
-        self, spatial_data: geopandas.GeoDataFrame
+        self, spatial_data: geopandas.GeoDataFrame, map_data: Optional[MapData] = None
     ) -> pandas.DataFrame:
         """
         Execute sample method takes full point data, samples the data and returns the decimated points
@@ -93,20 +87,10 @@ def sample(
         data = spatial_data.copy()
         data["X"] = data.geometry.x
         data["Y"] = data.geometry.y
-        if self.dtm_data is not None:
-            result = set_z_values_from_raster_df(self.dtm_data, data)
-            if result is not None:
-                data["Z"] = result["Z"]
-            else:
-                data["Z"] = None
-        else:
-            data["Z"] = None
-        if self.geology_data is not None:
-            data["layerID"] = geopandas.sjoin(
-                data, self.geology_data, how='left'
-            )['index_right']
-        else:
-            data["layerID"] = None
+        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']
         data.reset_index(drop=True, inplace=True)
 
         return pandas.DataFrame(data[:: self.decimation].drop(columns="geometry"))
@@ -121,24 +105,20 @@ class SamplerSpacing(Sampler):
     """
 
     @beartype.beartype
-    def __init__(self, spacing: float = 50.0, dtm_data: Optional[gdal.Dataset] = None, geology_data: Optional[geopandas.GeoDataFrame] = None):
+    def __init__(self, spacing: float = 50.0):
         """
         Initialiser for spacing sampler
 
         Args:
             spacing (float, optional): The distance between samples. Defaults to 50.0.
-            dtm_data (Optional[gdal.Dataset], optional): digital terrain map data. Defaults to None.
-            geology_data (Optional[geopandas.GeoDataFrame], optional): geology data. Defaults to None.
         """
-        super().__init__(dtm_data, geology_data)
         self.sampler_label = "SamplerSpacing"
         spacing = max(spacing, 1.0)
         self.spacing = spacing
 
-
     @beartype.beartype
     def sample(
-        self, spatial_data: geopandas.GeoDataFrame
+        self, spatial_data: geopandas.GeoDataFrame, map_data: Optional[MapData] = None
     ) -> pandas.DataFrame:
         """
         Execute sample method takes full point data, samples the data and returns the sampled points
diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py
index 4d25358c..d7a9aad1 100644
--- a/map2loop/thickness_calculator.py
+++ b/map2loop/thickness_calculator.py
@@ -5,17 +5,16 @@
     calculate_endpoints,
     multiline_to_line,
     find_segment_strike_from_pt,
-    set_z_values_from_raster_df,
-    value_from_raster
 )
+from .m2l_enums import Datatype
 from .interpolators import DipDipDirectionInterpolator
+from .mapdata import MapData
 
 from .logging import getLogger
 logger = getLogger(__name__)  
 
 # external imports
 from abc import ABC, abstractmethod
-from typing import Optional
 import scipy.interpolate
 import beartype
 import numpy
@@ -24,7 +23,6 @@
 import geopandas
 import shapely
 import math
-from osgeo import gdal
 from shapely.errors import UnsupportedGEOSVersionError
 
 class ThicknessCalculator(ABC):
@@ -35,19 +33,12 @@ class ThicknessCalculator(ABC):
         ABC (ABC): Derived from Abstract Base Class
     """
 
-    def __init__(
-        self, 
-        dtm_data: Optional[gdal.Dataset] = None, 
-        bounding_box: Optional[dict] = None, 
-        max_line_length: Optional[float] = None
-        ):
+    def __init__(self, max_line_length: float = None):
         """
         Initialiser of for ThicknessCalculator
         """
         self.thickness_calculator_label = "ThicknessCalculatorBaseClass"
         self.max_line_length = max_line_length
-        self.dtm_data = dtm_data
-        self.bounding_box = bounding_box
 
     def type(self):
         """
@@ -66,8 +57,7 @@ def compute(
         stratigraphic_order: list,
         basal_contacts: geopandas.GeoDataFrame,
         structure_data: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame,
-        sampled_contacts: pandas.DataFrame,
+        map_data: MapData,
     ) -> pandas.DataFrame:
         """
         Execute thickness calculator method (abstract method)
@@ -124,8 +114,7 @@ def compute(
         stratigraphic_order: list,
         basal_contacts: geopandas.GeoDataFrame,
         structure_data: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame,
-        sampled_contacts: pandas.DataFrame,
+        map_data: MapData,
     ) -> pandas.DataFrame:
         """
         Execute thickness calculator method takes unit data, basal_contacts and stratigraphic order and attempts to estimate unit thickness.
@@ -154,7 +143,7 @@ def compute(
 
         basal_unit_list = basal_contacts["basal_unit"].to_list()
         sampled_basal_contacts = rebuild_sampled_basal_contacts(
-            basal_contacts=basal_contacts, sampled_contacts=sampled_contacts
+            basal_contacts=basal_contacts, sampled_contacts=map_data.sampled_contacts
         )
 
         if len(stratigraphic_order) < 3:
@@ -217,16 +206,11 @@ class InterpolatedStructure(ThicknessCalculator):
         -> pandas.DataFrame: Calculates a thickness map for the overall map area.
     """
 
-    def __init__(
-        self, 
-        dtm_data: Optional[gdal.Dataset] = None, 
-        bounding_box: Optional[dict] = None, 
-        max_line_length: Optional[float] = None
-        ):
+    def __init__(self, max_line_length: float = None):
         """
         Initialiser for interpolated structure version of the thickness calculator
         """
-        super().__init__(dtm_data, bounding_box, max_line_length)
+        super().__init__(max_line_length)
         self.thickness_calculator_label = "InterpolatedStructure"
         self.lines = None
         
@@ -238,8 +222,7 @@ def compute(
         stratigraphic_order: list,
         basal_contacts: geopandas.GeoDataFrame,
         structure_data: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame,
-        sampled_contacts: pandas.DataFrame,
+        map_data: MapData,
     ) -> pandas.DataFrame:
         """
         Execute thickness calculator method takes unit data, basal_contacts, stratigraphic order, orientation data and
@@ -280,7 +263,7 @@ def compute(
         # increase buffer around basal contacts to ensure that the points are included as intersections
         basal_contacts["geometry"] = basal_contacts["geometry"].buffer(0.01)
         # get the sampled contacts
-        contacts = geopandas.GeoDataFrame(sampled_contacts)
+        contacts = geopandas.GeoDataFrame(map_data.sampled_contacts)
         # build points from x and y coordinates
         geometry2 = geopandas.points_from_xy(contacts['X'], contacts['Y'])
         contacts.set_geometry(geometry2, inplace=True)
@@ -288,7 +271,7 @@ 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 = set_z_values_from_raster_df(self.dtm_data, contacts)
+        contacts = map_data.get_value_from_raster_df(Datatype.DTM, contacts)
         # update the geometry of the contact points to include the Z value
         contacts["geometry"] = contacts.apply(
             lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
@@ -296,10 +279,12 @@ def compute(
         # spatial join the contact points with the basal contacts to get the unit for each contact point
         contacts = contacts.sjoin(basal_contacts, how="inner", predicate="intersects")
         contacts = contacts[["X", "Y", "Z", "geometry", "basal_unit"]].copy()
+        bounding_box = map_data.get_bounding_box()
+
         # Interpolate the dip of the contacts
         interpolator = DipDipDirectionInterpolator(data_type="dip")
         # Interpolate the dip of the contacts
-        dip = interpolator(self.bounding_box, structure_data, interpolator=scipy.interpolate.Rbf)
+        dip = interpolator(bounding_box, structure_data, interpolator=scipy.interpolate.Rbf)
         # create a GeoDataFrame of the interpolated orientations
         interpolated_orientations = geopandas.GeoDataFrame()
         # add the dip and dip direction to the GeoDataFrame
@@ -314,13 +299,13 @@ 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 = set_z_values_from_raster_df(self.dtm_data, interpolated_orientations)
+        interpolated = map_data.get_value_from_raster_df(Datatype.DTM, interpolated_orientations)
         # update the geometry of the interpolated points to include the Z value
         interpolated["geometry"] = interpolated.apply(
             lambda row: shapely.geometry.Point(row.geometry.x, row.geometry.y, row["Z"]), axis=1
         )
         # for each interpolated point, assign name of unit using spatial join
-        units = geology_data.copy()
+        units = map_data.get_map_data(Datatype.GEOLOGY)
         interpolated_orientations = interpolated_orientations.sjoin(
             units, how="inner", predicate="within"
         )
@@ -364,22 +349,19 @@ def compute(
                         # check if the short line is 
                         if self.max_line_length is not None and short_line.length > self.max_line_length:
                             continue
-                        
-                        inv_geotransform = gdal.InvGeoTransform(self.dtm_data.GetGeoTransform())
-                        data_array = numpy.array(self.dtm_data.GetRasterBand(1).ReadAsArray().T)
  
                         # extract the end points of the shortest line
                         p1 = numpy.zeros(3)
                         p1[0] = numpy.asarray(short_line[0].coords[0][0])
                         p1[1] = numpy.asarray(short_line[0].coords[0][1])
                         # get the elevation Z of the end point p1
-                        p1[2] = value_from_raster(inv_geotransform, data_array, p1[0], p1[1])
+                        p1[2] = map_data.get_value_from_raster(Datatype.DTM, p1[0], p1[1])
                         # create array to store xyz coordinates of the end point p2
                         p2 = numpy.zeros(3)
                         p2[0] = numpy.asarray(short_line[0].coords[-1][0])
                         p2[1] = numpy.asarray(short_line[0].coords[-1][1])
                         # get the elevation Z of the end point p2
-                        p2[2] = value_from_raster(inv_geotransform, data_array, p2[0], p2[1])
+                        p2[2] = map_data.get_value_from_raster(Datatype.DTM, p2[0], p2[1])
                         # calculate the length of the shortest line
                         line_length = scipy.spatial.distance.euclidean(p1, p2)
                         # find the indices of the points that are within 5% of the length of the shortest line
@@ -459,13 +441,8 @@ class StructuralPoint(ThicknessCalculator):
 
     '''
 
-    def __init__(
-        self, 
-        dtm_data: Optional[gdal.Dataset] = None, 
-        bounding_box: Optional[dict] = None, 
-        max_line_length: Optional[float] = None
-        ):
-        super().__init__(dtm_data, bounding_box, max_line_length)
+    def __init__(self, max_line_length: float = None):
+        super().__init__(max_line_length)
         self.thickness_calculator_label = "StructuralPoint"
         self.strike_allowance = 30
         self.lines = None
@@ -478,8 +455,7 @@ def compute(
         stratigraphic_order: list,
         basal_contacts: geopandas.GeoDataFrame,
         structure_data: pandas.DataFrame,
-        geology_data: geopandas.GeoDataFrame,
-        sampled_contacts: pandas.DataFrame,
+        map_data: MapData,
     ) -> pandas.DataFrame:
         """
         Method overview:
@@ -520,7 +496,7 @@ def compute(
         basal_contacts = basal_contacts.copy()
 
         # grab geology polygons and calculate bounding boxes for each lithology
-        geology = geology_data.copy()
+        geology = map_data.get_map_data(datatype=Datatype.GEOLOGY)
         geology[['minx', 'miny', 'maxx', 'maxy']] = geology.bounds
 
         # create a GeoDataFrame of the sampled structures
@@ -541,7 +517,7 @@ def compute(
 
         # rebuild basal contacts lines based on sampled dataset
         sampled_basal_contacts = rebuild_sampled_basal_contacts(
-            basal_contacts, sampled_contacts
+            basal_contacts, map_data.sampled_contacts
         )
 
         # calculate map dimensions
diff --git a/map2loop/utils.py b/map2loop/utils.py
index 55e2e7b2..c3ed7795 100644
--- a/map2loop/utils.py
+++ b/map2loop/utils.py
@@ -7,7 +7,6 @@
 import pandas
 import re
 import json
-from osgeo import gdal
 
 from .logging import getLogger
 logger = getLogger(__name__)
@@ -529,62 +528,3 @@ def update_from_legacy_file(
             json.dump(parsed_data, f, indent=4)
     
     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
\ No newline at end of file
diff --git a/tests/contact_extractor/test_contact_extractor.py b/tests/contact_extractor/test_contact_extractor.py
deleted file mode 100644
index 1c5acef6..00000000
--- a/tests/contact_extractor/test_contact_extractor.py
+++ /dev/null
@@ -1,37 +0,0 @@
-import sys
-sys.path.append('/usr/lib/python3/dist-packages')
-from map2loop.contact_extractor import ContactExtractor
-import geopandas as gpd
-from shapely.geometry import Polygon
-
-def simple_geology():
-    poly1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
-    poly2 = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)])
-    return gpd.GeoDataFrame(
-        {
-            "UNITNAME": ["A", "B"],
-            "INTRUSIVE": [False, False],
-            "SILL": [False, False],
-            "geometry": [poly1, poly2],
-        },
-        geometry="geometry",
-        crs="EPSG:28350",
-    )
-
-def test_extract_all_contacts():
-    geology = simple_geology()
-    extractor = ContactExtractor(geology, None)
-    contacts = extractor.extract_all_contacts()
-    assert len(contacts) == 1
-    assert {contacts.loc[0, "UNITNAME_1"], contacts.loc[0, "UNITNAME_2"]} == {"A", "B"}
-
-def test_extract_basal_contacts():
-    geology = simple_geology()
-    extractor = ContactExtractor(geology, None)
-    extractor.extract_all_contacts()
-    basal = extractor.extract_basal_contacts(["A", "B"], save_contacts=True)
-    assert len(basal) == 1
-    assert basal.loc[0, "basal_unit"] == "A"
-    assert basal.loc[0, "type"] == "BASAL"
-    assert extractor.basal_contacts is not None
-    assert extractor.all_basal_contacts is not None
diff --git a/tests/project/test_plot_hamersley.py b/tests/project/test_plot_hamersley.py
index 504c4585..07393f27 100644
--- a/tests/project/test_plot_hamersley.py
+++ b/tests/project/test_plot_hamersley.py
@@ -31,13 +31,13 @@ def create_project(state_data="WA", projection="EPSG:28350"):
 
 # is the project running?
 def test_project_execution():
-    try:
-        proj = create_project()
-    except Exception:
-        pytest.skip("Skipping the project test from server data due to loading failure")
+
+    proj = create_project()
     try:
         proj.run_all(take_best=True)
+    # if there's a timeout:
     except requests.exceptions.ReadTimeout:
+        print("Timeout occurred, skipping the test.")  # Debugging line
         pytest.skip(
             "Skipping the project test from server data due to timeout while attempting to run proj.run_all"
         )
diff --git a/tests/project/test_thickness_calculations.py b/tests/project/test_thickness_calculations.py
index 373687ed..9cf5464d 100644
--- a/tests/project/test_thickness_calculations.py
+++ b/tests/project/test_thickness_calculations.py
@@ -1706,14 +1706,14 @@ def test_calculate_unit_thicknesses():
 
     # check set
 
-    project.set_thickness_calculator([StructuralPoint(dtm_data=dtm, bounding_box=bbox_3d), InterpolatedStructure(dtm_data=dtm, bounding_box=bbox_3d)])
+    project.set_thickness_calculator([StructuralPoint(), InterpolatedStructure()])
     assert project.get_thickness_calculator() == [
         'StructuralPoint',
         'InterpolatedStructure',
     ], "Setter method for thickness calculator not working"  ## default is InterpolatedStructure
 
     # Run the calculate_unit_thicknesses
-    project.calculate_unit_thicknesses(bc_gdf)
+    project.calculate_unit_thicknesses()
 
     # # Check if all thicknesses have been calculated
     columns_to_check = [
diff --git a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py
index 12d0c87c..9bc28431 100644
--- a/tests/thickness/InterpolatedStructure/test_interpolated_structure.py
+++ b/tests/thickness/InterpolatedStructure/test_interpolated_structure.py
@@ -1663,6 +1663,7 @@ def check_thickness_values(result, column, description):
 
 def test_calculate_thickness_InterpolatedStructure():
     # Run the calculation
+    thickness_calculator = InterpolatedStructure()
 
     md = MapData()
     md.sampled_contacts = s_c
@@ -1682,16 +1683,13 @@ def test_calculate_thickness_InterpolatedStructure():
         "base": -3200,
         "top": 3000,
     }
-    thickness_calculator = InterpolatedStructure(dtm_data=md.get_map_data(Datatype.DTM),
-                                                 bounding_box=md.bounding_box)
 
     result = thickness_calculator.compute(
         units=st_units,
         stratigraphic_order=st_column,
         basal_contacts=bc_gdf,
         structure_data=structures,
-        geology_data=md.get_map_data(Datatype.GEOLOGY),
-        sampled_contacts=md.sampled_contacts,
+        map_data=md,
     )
 
     # is thickness calc alpha the label?
diff --git a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py
index f528b0ee..a71ebb7b 100644
--- a/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py
+++ b/tests/thickness/StructurePoint/test_ThicknessStructuralPoint.py
@@ -1664,6 +1664,7 @@ def test_calculate_thickness_structural_point():
 
     md = MapData()
     md.sampled_contacts = s_c
+    md.sampled_contacts = s_c
     md.raw_data[Datatype.GEOLOGY] = geology
     md.load_map_data(Datatype.GEOLOGY)
     md.check_map(Datatype.GEOLOGY)
@@ -1674,8 +1675,7 @@ def test_calculate_thickness_structural_point():
         stratigraphic_order=st_column,
         basal_contacts=bc_gdf,
         structure_data=structures,
-        geology_data=md.get_map_data(Datatype.GEOLOGY),
-        sampled_contacts=md.sampled_contacts,
+        map_data=md,
     )
 
     # is thickness calc alpha the label?
diff --git a/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py b/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py
index 521b633f..a169b3ce 100644
--- a/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py
+++ b/tests/thickness/ThicknessCalculatorAlpha/test_ThicknessCalculatorAlpha.py
@@ -1641,6 +1641,7 @@
 
 s_c = pandas.DataFrame({'X': X, 'Y': Y, 'Z': Z, 'featureId': featureid})
 
+
 #####################################
 ### TEST ThicknessCalculatorAlpha ###
 #####################################
@@ -1657,16 +1658,15 @@ def check_thickness_values(result, column, description):
 
 
 geology = load_hamersley_geology()
-geology.rename(columns={'unitname': 'UNITNAME', 'code': 'CODE'}, inplace=True)
+
 
 def test_calculate_thickness_thickness_calculator_alpha():
     # Run the calculation
     md = MapData()
     md.sampled_contacts = s_c
-    md.raw_data[Datatype.GEOLOGY] = geology
-    md.load_map_data(Datatype.GEOLOGY)
-    md.check_map(Datatype.GEOLOGY)
-    md.parse_geology_map()
+    md.data[Datatype.GEOLOGY] = geology
+
+    print('GERE', md.get_map_data(Datatype.GEOLOGY))
 
     thickness_calculator = ThicknessCalculatorAlpha()
     result = thickness_calculator.compute(
@@ -1674,8 +1674,7 @@ def test_calculate_thickness_thickness_calculator_alpha():
         stratigraphic_order=st_column,
         basal_contacts=bc_gdf,
         structure_data=structures,
-        geology_data=md.get_map_data(Datatype.GEOLOGY),
-        sampled_contacts=md.sampled_contacts,
+        map_data=md,
     )
 
     # is thickness calc alpha the label?
diff --git a/tests/topology/test_topology.py b/tests/topology/test_topology.py
deleted file mode 100644
index d5adce10..00000000
--- a/tests/topology/test_topology.py
+++ /dev/null
@@ -1,90 +0,0 @@
-import geopandas as gpd
-from shapely.geometry import Polygon, LineString
-
-from map2loop.topology import Topology
-
-
-
-def simple_geology():
-    poly1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
-    poly2 = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)])
-    return gpd.GeoDataFrame(
-        {
-            "UNITNAME": ["A", "B"],
-            "INTRUSIVE": [False, False],
-            "SILL": [False, False],
-            "geometry": [poly1, poly2],
-        },
-        geometry="geometry",
-        crs="EPSG:28350",
-    )
-
-
-def faults():
-    line1 = LineString([(0, 0), (2, 1)])
-    line2 = LineString([(0, 1), (2, 0)])
-    return gpd.GeoDataFrame(
-        {"ID": [1, 2], "geometry": [line1, line2]}, geometry="geometry", crs="EPSG:28350"
-    )
-
-
-def test_initialisation_defaults():
-    geology_data = simple_geology()
-    faults_data = faults()
-    topo = Topology(geology_data, faults_data)
-
-    assert topo.sorted_units is None
-    assert topo._fault_fault_relationships is None
-    assert topo._unit_fault_relationships is None
-    assert topo._unit_unit_relationships is None
-    assert topo.buffer_radius == 500
-
-
-def test_calculate_fault_fault_relationships():
-    geology_data = simple_geology()
-    faults_data = faults()
-    topo = Topology(geology_data, faults_data)
-    topo.buffer_radius = 0.1
-
-    df = topo.get_fault_fault_relationships()
-    assert list(df.columns) == ["Fault1", "Fault2", "Angle", "Type"]
-    assert len(df) == 1
-    assert set(df.loc[0, ["Fault1", "Fault2"]]) == {1, 2}
-    assert df.loc[0, "Angle"] == 60
-    assert df.loc[0, "Type"] == "T"
-
-
-def test_calculate_unit_fault_relationships():
-    geology_data = simple_geology()
-    faults_data = faults()
-    topo = Topology(geology_data, faults_data)
-    topo.buffer_radius = 0.1
-
-    df = topo.get_unit_fault_relationships()
-    assert set(df["Unit"]) == {"A", "B"}
-    assert set(df["Fault"]) == {1, 2}
-    # each unit is intersected by both faults
-    assert len(df) == 4
-
-
-def test_calculate_unit_unit_relationships_with_contacts():
-    geology_data = simple_geology()
-    topo = Topology(geology_data, None)
-    df = topo.get_unit_unit_relationships()
-    assert list(df.columns) == ["UNITNAME_1", "UNITNAME_2"]
-    assert len(df) == 1
-    assert set(df.loc[0]) == {"A", "B"}
-
-
-def test_reset():
-    geology_data = simple_geology()
-    faults_data = faults()
-    topo = Topology(geology_data, faults_data)
-    topo.buffer_radius = 0.1
-    ffr = topo.get_fault_fault_relationships()
-    assert ffr is not None
-
-    topo.reset()
-    assert topo._fault_fault_relationships is None
-    assert topo._unit_fault_relationships is None
-    assert topo._unit_unit_relationships is None
\ No newline at end of file