diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 36845786..4cf1d570 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -31,11 +31,11 @@ jobs:
       - name: checkout code
         uses: actions/checkout@v3
 
-      - name: Setup Mambaforge
+      - name: Setup Miniforge
         uses: conda-incubator/setup-miniconda@v2
         with:
             python-version: '3.11'
-            miniforge-variant: Mambaforge
+            miniforge-variant: Miniforge3
             miniforge-version: latest
             use-mamba: true
 
diff --git a/.github/workflows/tests_dev.yml b/.github/workflows/tests_dev.yml
index d67ae839..fcb12dac 100644
--- a/.github/workflows/tests_dev.yml
+++ b/.github/workflows/tests_dev.yml
@@ -1,62 +1,63 @@
-name: Tests with HydroMT dev
+# Uncomment when we move to HydroMT v1
+# name: Tests with HydroMT dev
 
-on:
-  # trigger weekly on monday at 00:00 UTC
-  schedule:
-    - cron: '0 0 * * 1'
-  push:
-    branches: [main]
-    paths:
-      - .github/workflows/tests_dev.yml
-      - tests/*
-      - hydromt_sfincs/*
-      - pyproject.toml
-  pull_request:
-    branches: [main]
-    paths:
-      - .github/workflows/tests_dev.yml
-      - tests/*
-      - hydromt_sfincs/*
-      - pyproject.toml
+# on:
+#   # trigger weekly on monday at 00:00 UTC
+#   schedule:
+#     - cron: '0 0 * * 1'
+#   push:
+#     branches: [main]
+#     paths:
+#       - .github/workflows/tests_dev.yml
+#       - tests/*
+#       - hydromt_sfincs/*
+#       - pyproject.toml
+#   pull_request:
+#     branches: [main]
+#     paths:
+#       - .github/workflows/tests_dev.yml
+#       - tests/*
+#       - hydromt_sfincs/*
+#       - pyproject.toml
 
-jobs:
-  build:
-    defaults:
-      run:
-        shell: bash -l {0}
-    strategy:
-      fail-fast: false
-    runs-on: ubuntu-latest
-    timeout-minutes: 30
-    concurrency:
-      group: ${{ github.workflow }}-${{ matrix.python-version }}-${{ github.ref }}
-      cancel-in-progress: true
-    steps:
+# jobs:
+#   build:
+#     defaults:
+#       run:
+#         shell: bash -l {0}
+#     strategy:
+#       fail-fast: false
+#     runs-on: ubuntu-latest
+#     timeout-minutes: 30
+#     concurrency:
+#       group: ${{ github.workflow }}-${{ matrix.python-version }}-${{ github.ref }}
+#       cancel-in-progress: true
+#     steps:
 
-      - uses: actions/checkout@v3
+#       - uses: actions/checkout@v3
 
-      - uses: actions/setup-python@v5
-        id: pip
-        with:
-          # caching, see https://github.com/actions/setup-python/blob/main/docs/advanced-usage.md#caching-packages
-          cache: 'pip'
-          python-version: '3.9'  # 3.9 is not supported by last released version of hydromt
-          cache-dependency-path: pyproject.toml
+#       - uses: actions/setup-python@v5
+#         id: pip
+#         with:
+#           # caching, see https://github.com/actions/setup-python/blob/main/docs/advanced-usage.md#caching-packages
+#           cache: 'pip'
+#           python-version: '3.9'  # 3.9 is not supported by last released version of hydromt
+#           cache-dependency-path: pyproject.toml
 
-      # true if cache-hit occurred on the primary key
-      - name: Cache hit
-        run: echo '${{ steps.pip.outputs.cache-hit }}'
+#       # true if cache-hit occurred on the primary key
+#       - name: Cache hit
+#         run: echo '${{ steps.pip.outputs.cache-hit }}'
 
-      # build environment with pip
-      - name: Install hydromt-sfincs
-        run: |
-          pip install --upgrade pip
-          pip install .[test,examples]
-          pip install git+https://github.com/Deltares/hydromt.git
-          pip list
+#       # build environment with pip
+#       - name: Install hydromt-sfincs
+#         run: |
+#           pip install --upgrade pip
+#           pip install .[test,examples]
+#           pip install git+https://github.com/Deltares/hydromt.git
+#           pip list
 
-      # run test
-      - name: Test
-        run: |
-          export NUMBA_DISABLE_JIT=1
-          python -m pytest --verbose
+#       # run test
+#       - name: Test
+#         run: |
+#           export NUMBA_DISABLE_JIT=1
+#           python -m pytest --verbose
diff --git a/envs/hydromt-sfincs-dev.yml b/envs/hydromt-sfincs-dev.yml
index aede86bf..6a23c867 100644
--- a/envs/hydromt-sfincs-dev.yml
+++ b/envs/hydromt-sfincs-dev.yml
@@ -6,6 +6,7 @@ channels:
 dependencies:
 - affine
 - cartopy
+- datashader
 - descartes
 - ffmpeg
 - geopandas>1.0
@@ -30,6 +31,7 @@ dependencies:
 - shapely
 - sphinx
 - xarray
+- xugrid
 - pip:
   - black[jupyter]
   - sphinx_design
diff --git a/hydromt_sfincs/plots.py b/hydromt_sfincs/plots.py
index 7859ae8e..e5be2c47 100644
--- a/hydromt_sfincs/plots.py
+++ b/hydromt_sfincs/plots.py
@@ -1,13 +1,16 @@
 """Plotting functions for SFINCS model data."""
-
-from typing import Dict, List, Tuple
+from typing import Dict, List, Tuple, Union
 
 import numpy as np
+import logging
 import pandas as pd
 import xarray as xr
+import xugrid as xu
 
 from .utils import get_bounds_vector
 
+logger = logging.getLogger(__name__)
+
 __all__ = ["plot_forcing", "plot_basemap"]
 
 geom_style = {
@@ -107,7 +110,7 @@ def plot_forcing(forcing: Dict, **kwargs):
 
 
 def plot_basemap(
-    ds: xr.Dataset,
+    ds: Union[xr.Dataset, xu.UgridDataset],
     geoms: Dict,
     variable: str = "dep",
     shaded: bool = False,
@@ -121,6 +124,7 @@ def plot_basemap(
     geom_kwargs: Dict = {},
     legend_kwargs: Dict = {},
     bmap_kwargs: Dict = {},
+    logger=logger,
     **kwargs,
 ):
     """Create basemap plot.
@@ -175,10 +179,43 @@ def plot_basemap(
         has_cx = False
 
     # read crs and utm zone > convert to cartopy
-    proj_crs = ds.raster.crs
+    if isinstance(ds, xr.Dataset):
+        proj_crs = ds.raster.crs
+        bounds = ds.raster.box.buffer(abs(ds.raster.res[0] * 1)).total_bounds
+        bbox = ds.raster.transform_bounds(4326)
+        ratio = ds.raster.ycoords.size / (ds.raster.xcoords.size * 1.4)
+    elif isinstance(ds, xu.UgridDataset):
+        proj_crs = ds.grid.crs
+        bounds = ds.ugrid.total_bounds
+        bbox = ds.ugrid.to_crs(4326).ugrid.total_bounds
+        ratio = (bounds[3] - bounds[1]) / ((bounds[2] - bounds[0]) * 1.4)
     proj_str = proj_crs.name
+    extent = np.array(bounds)[[0, 2, 1, 3]]
+
     if proj_crs.is_projected and proj_crs.to_epsg() is not None:
-        crs = ccrs.epsg(ds.raster.crs.to_epsg())
+        if "UTM" in proj_str:
+            # Extract the UTM zone number
+            utm_zone = proj_crs.utm_zone
+            # Parse the zone number and hemisphere
+            zone_number = int(utm_zone[:-1])  # Extract numeric part
+            hemisphere = utm_zone[-1]  # Last character is 'N' or 'S'
+            # Determine if it's in the southern hemisphere
+            is_southern_hemisphere = hemisphere == "S"
+            # Create the Cartopy UTM projection
+            crs = ccrs.UTM(zone_number, southern_hemisphere=is_southern_hemisphere)
+        else:
+            crs = ccrs.epsg(proj_crs.to_epsg())
+
+        # now check whether the model extent is within the extent of the crs
+        bounds_crs = crs.boundary.bounds
+        if (
+            bounds[0] < bounds_crs[0]
+            or bounds[1] < bounds_crs[1]
+            or bounds[2] > bounds_crs[2]
+            or bounds[3] > bounds_crs[3]
+        ):
+            logger.warning("Model domain exceeds the area of the CRS.")
+            logger.warning("Consider reprojecting to EPSG:4326 for plotting.")
         unit = proj_crs.axis_info[0].unit_name
         unit = "m" if unit == "metre" else unit
         xlab, ylab = f"x [{unit}] - {proj_str}", f"y [{unit}]"
@@ -187,12 +224,9 @@ def plot_basemap(
         xlab, ylab = f"lon [deg] - {proj_str}", "lat [deg]"
     else:
         raise ValueError("Unsupported CRS")
-    bounds = ds.raster.box.buffer(abs(ds.raster.res[0] * 1)).total_bounds
-    extent = np.array(bounds)[[0, 2, 1, 3]]
 
     # create fig with geo-axis and set background
     if figsize is None:
-        ratio = ds.raster.ycoords.size / (ds.raster.xcoords.size * 1.4)
         figsize = (6, 6 * ratio)
     fig = plt.figure(figsize=figsize)
     ax = plt.subplot(projection=crs)
@@ -200,7 +234,7 @@ def plot_basemap(
     if bmap is not None:
         if zoomlevel == "auto":  # auto zoomlevel
             c = 2 * np.pi * 6378137  # Earth circumference
-            lat = np.array(ds.raster.transform_bounds(4326))[[1, 3]].mean()
+            lat = np.array(bbox)[[1, 3]].mean()
             # max 4 x 4 tiles per image
             tile_size = max(bounds[2] - bounds[0], bounds[3] - bounds[1]) / 4
             if proj_crs.is_geographic:
@@ -247,30 +281,35 @@ def plot_basemap(
             kwargs0["cbar_kwargs"].update(ticks=[1, 2, 3])
 
     if variable in ds:
-        da = ds[variable].raster.mask_nodata()
+        da = ds[variable]
         if "msk" in ds and np.any(ds["msk"] > 0):
             da = da.where(ds["msk"] > 0)
-        if da.raster.rotation != 0 and "xc" in da.coords and "yc" in da.coords:
-            da.plot(transform=crs, x="xc", y="yc", ax=ax, zorder=1, **kwargs0)
-        else:
-            da.plot.imshow(transform=crs, ax=ax, zorder=1, **kwargs0)
-        if shaded and variable == "dep" and da.raster.rotation == 0:
-            ls = colors.LightSource(azdeg=315, altdeg=45)
-            dx, dy = da.raster.res
-            _rgb = ls.shade(
-                da.fillna(0).values,
-                norm=kwargs0["norm"],
-                cmap=kwargs0["cmap"],
-                blend_mode="soft",
-                dx=dx,
-                dy=dy,
-                vert_exag=2,
-            )
-            rgb = xr.DataArray(
-                dims=("y", "x", "rgb"), data=_rgb, coords=da.raster.coords
-            )
-            rgb = xr.where(np.isnan(da), np.nan, rgb)
-            rgb.plot.imshow(transform=crs, ax=ax, zorder=1)
+        if isinstance(da, xr.DataArray):
+            # mask_nodata converts it to an xr.DataArray, so performed inside if-else statement
+            da = da.raster.mask_nodata()
+            if da.raster.rotation != 0 and "xc" in da.coords and "yc" in da.coords:
+                da.plot(transform=crs, x="xc", y="yc", ax=ax, zorder=1, **kwargs0)
+            else:
+                da.plot.imshow(transform=crs, ax=ax, zorder=1, **kwargs0)
+            if shaded and variable == "dep" and da.raster.rotation == 0:
+                ls = colors.LightSource(azdeg=315, altdeg=45)
+                dx, dy = da.raster.res
+                _rgb = ls.shade(
+                    da.fillna(0).values,
+                    norm=kwargs0["norm"],
+                    cmap=kwargs0["cmap"],
+                    blend_mode="soft",
+                    dx=dx,
+                    dy=dy,
+                    vert_exag=2,
+                )
+                rgb = xr.DataArray(
+                    dims=("y", "x", "rgb"), data=_rgb, coords=da.raster.coords
+                )
+                rgb = xr.where(np.isnan(da), np.nan, rgb)
+                rgb.plot.imshow(transform=crs, ax=ax, zorder=1)
+        elif isinstance(da, xu.UgridDataArray):
+            da.ugrid.plot(transform=crs, ax=ax, zorder=1, **kwargs0)
 
     # geometry plotting and annotate kwargs
     for k, d in geom_kwargs.items():
@@ -290,6 +329,11 @@ def plot_basemap(
             "No 'msk' (sfincs.msk) found in ds required to plot the model bounds "
             "Set plot_bounds=False or add 'msk' to ds"
         )
+    elif plot_bounds and isinstance(ds, xu.UgridDataset):
+        raise NotImplementedError(
+            "Plotting of the boundaries for quadtree grids is not yet implemented. "
+            "Set plot_bounds=False to proceed."
+        )
     elif plot_bounds and (ds["msk"] >= 1).any():
         gdf_msk = get_bounds_vector(ds["msk"])
         gdf_msk2 = gdf_msk[gdf_msk["value"] == 2]
diff --git a/hydromt_sfincs/quadtree.py b/hydromt_sfincs/quadtree.py
new file mode 100644
index 00000000..94086ba5
--- /dev/null
+++ b/hydromt_sfincs/quadtree.py
@@ -0,0 +1,192 @@
+import logging
+import os
+from pathlib import Path
+from typing import Union
+
+import numpy as np
+from pyproj import CRS, Transformer
+
+import geopandas as gpd
+import pandas as pd
+import xarray as xr
+import xugrid as xu
+import shapely
+
+# optional dependency
+try:
+    from datashader import Canvas
+    import datashader.transfer_functions as tf
+    from datashader.utils import export_image
+
+    HAS_DATASHADER = True
+except ImportError:
+    HAS_DATASHADER = False
+
+
+from hydromt_sfincs.subgrid import SubgridTableQuadtree
+
+logger = logging.getLogger(__name__)
+
+
+class QuadtreeGrid:
+    def __init__(self, logger=logger):
+        self.nr_cells = 0
+        self.nr_refinement_levels = 1
+        self.version = 0
+
+        self.data = None  # placeholder for xugrid object
+        self.subgrid = SubgridTableQuadtree()
+        self.df = None  # placeholder for pandas dataframe for datashader
+
+    @property
+    def crs(self):
+        if self.data is None:
+            return None
+        return self.data.grid.crs
+
+    @property
+    def face_coordinates(self):
+        if self.data is None:
+            return None
+        xy = self.data.grid.face_coordinates
+        return xy[:, 0], xy[:, 1]
+
+    @property
+    def exterior(self):
+        if self.data is None:
+            return gpd.GeoDataFrame()
+        indx = self.data.grid.edge_node_connectivity[self.data.grid.exterior_edges, :]
+        x = self.data.grid.node_x[indx]
+        y = self.data.grid.node_y[indx]
+
+        # Make linestrings from numpy arrays x and y
+        linestrings = [
+            shapely.LineString(np.column_stack((x[i], y[i]))) for i in range(len(x))
+        ]
+        # Merge linestrings
+        merged = shapely.ops.linemerge(linestrings)
+        # Merge polygons
+        polygons = shapely.ops.polygonize(merged)
+
+        return gpd.GeoDataFrame(geometry=list(polygons), crs=self.crs)
+
+    @property
+    def empty_mask(self):
+        if self.data is None:
+            return None
+        # create empty mask
+        da0 = xr.DataArray(
+            data=np.zeros(shape=len(self.data.grid.face_coordinates)),
+            dims=self.data.grid.face_dimension,
+        )
+        return xu.UgridDataArray(da0, self.data.grid)
+
+    def read(self, file_name: Union[str, Path] = "sfincs.nc"):
+        """Reads a quadtree netcdf file and stores it in the QuadtreeGrid object."""
+
+        self.data = xu.open_dataset(file_name)
+        self.data.close()  # TODO check if close works/is needed
+
+        # TODO make similar to fortran conventions?
+        # Rename to python conventions
+        self.data = self.data.rename({"z": "dep"}) if "z" in self.data else self.data
+        self.data = (
+            self.data.rename({"mask": "msk"}) if "mask" in self.data else self.data
+        )
+        self.data = (
+            self.data.rename({"snapwave_mask": "snapwave_msk"})
+            if "snapwave_mask" in self.data
+            else self.data
+        )
+
+        self.nr_cells = self.data.sizes["mesh2d_nFaces"]
+
+        # set CRS (not sure if that should be stored in the netcdf in this way)
+        # self.data.crs = CRS.from_wkt(self.data["crs"].crs_wkt)
+        self.data.grid.set_crs(CRS.from_wkt(self.data["crs"].crs_wkt))
+
+        for key, value in self.data.attrs.items():
+            setattr(self, key, value)
+
+    def write(self, file_name: Union[str, Path] = "sfincs.nc", version: int = 0):
+        """Writes a quadtree SFINCS netcdf file."""
+
+        # TODO do we want to cut inactive cells here? Or already when creating the mask?
+
+        attrs = self.data.attrs
+        ds = self.data.ugrid.to_dataset()
+
+        # TODO make similar to fortran conventions
+        # RENAME TO FORTRAN CONVENTION
+        ds = ds.rename({"dep": "z"}) if "dep" in ds else ds
+        ds = ds.rename({"msk": "mask"}) if "msk" in ds else ds
+        ds = (
+            ds.rename({"snapwave_msk": "snapwave_mask"}) if "snapwave_msk" in ds else ds
+        )
+
+        ds.attrs = attrs
+        ds.to_netcdf(file_name)
+
+    def map_overlay(self, file_name, xlim=None, ylim=None, color="black", width=800):
+        # check if datashader is available
+        if not HAS_DATASHADER:
+            logger.warning("Datashader is not available. Please install datashader.")
+            return False
+        if self.data is None:
+            # No grid (yet)
+            return False
+        try:
+            if not hasattr(self, "df"):
+                self.df = None
+            if self.df is None:
+                self._get_datashader_dataframe()
+
+            transformer = Transformer.from_crs(4326, 3857, always_xy=True)
+            xl0, yl0 = transformer.transform(xlim[0], ylim[0])
+            xl1, yl1 = transformer.transform(xlim[1], ylim[1])
+            xlim = [xl0, xl1]
+            ylim = [yl0, yl1]
+            ratio = (ylim[1] - ylim[0]) / (xlim[1] - xlim[0])
+            height = int(width * ratio)
+            cvs = Canvas(
+                x_range=xlim, y_range=ylim, plot_height=height, plot_width=width
+            )
+            agg = cvs.line(self.df, x=["x1", "x2"], y=["y1", "y2"], axis=1)
+            img = tf.shade(agg)
+            path = os.path.dirname(file_name)
+            if not path:
+                path = os.getcwd()
+            name = os.path.basename(file_name)
+            name = os.path.splitext(name)[0]
+            export_image(img, name, export_path=path)
+            return True
+        except Exception as e:
+            logger.warning("Failed to create map overlay. Error: %s" % e)
+            return False
+
+    def snap_to_grid(self, polyline, max_snap_distance=1.0):
+        if len(polyline) == 0:
+            return gpd.GeoDataFrame()
+        geom_list = []
+        for _, line in polyline.iterrows():
+            geom = line["geometry"]
+            if geom.geom_type == "LineString":
+                geom_list.append(geom)
+        gdf = gpd.GeoDataFrame({"geometry": geom_list})
+        _, snapped_gdf = xu.snap_to_grid(
+            gdf, self.data.grid, max_snap_distance=max_snap_distance
+        )
+        snapped_gdf = snapped_gdf.set_crs(self.crs)
+        return snapped_gdf
+
+    # Internal functions
+    def _get_datashader_dataframe(self):
+        # Create a dataframe with line elements
+        x1 = self.data.grid.edge_node_coordinates[:, 0, 0]
+        x2 = self.data.grid.edge_node_coordinates[:, 1, 0]
+        y1 = self.data.grid.edge_node_coordinates[:, 0, 1]
+        y2 = self.data.grid.edge_node_coordinates[:, 1, 1]
+        transformer = Transformer.from_crs(self.crs, 3857, always_xy=True)
+        x1, y1 = transformer.transform(x1, y1)
+        x2, y2 = transformer.transform(x2, y2)
+        self.df = pd.DataFrame(dict(x1=x1, y1=y1, x2=x2, y2=y2))
diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py
index a1b24a89..590fc386 100644
--- a/hydromt_sfincs/sfincs.py
+++ b/hydromt_sfincs/sfincs.py
@@ -16,6 +16,8 @@
 import numpy as np
 import pandas as pd
 import xarray as xr
+import xugrid as xu
+from xugrid.core.wrap import UgridDataArray
 from hydromt.models.model_grid import GridModel
 from hydromt.vector import GeoDataArray, GeoDataset
 from hydromt.workflows.forcing import da_to_timedelta
@@ -24,6 +26,7 @@
 
 from . import DATADIR, plots, utils, workflows
 from .regulargrid import RegularGrid
+from .quadtree import QuadtreeGrid
 from .subgrid import SubgridTableRegular
 from .sfincs_input import SfincsInput
 
@@ -139,6 +142,16 @@ def __init__(
         self.quadtree = None
         self.subgrid = xr.Dataset()
 
+    def __del__(self):
+        """Close the model and remove the logger file handler."""
+        for handler in self.logger.handlers:
+            if (
+                isinstance(handler, logging.FileHandler)
+                and "hydromt.log" in handler.baseFilename
+            ):
+                handler.close()
+                self.logger.removeHandler(handler)
+
     @property
     def mask(self) -> xr.DataArray | None:
         """Returns model mask"""
@@ -147,6 +160,11 @@ def mask(self) -> xr.DataArray | None:
                 return self.grid["msk"]
             elif self.reggrid is not None:
                 return self.reggrid.empty_mask
+        elif self.grid_type == "quadtree":
+            if "msk" in self.quadtree.data:
+                return self.quadtree.data["msk"]
+            elif self.quadtree is not None:
+                return self.quadtree.empty_mask
 
     @property
     def region(self) -> gpd.GeoDataFrame:
@@ -155,21 +173,40 @@ def region(self) -> gpd.GeoDataFrame:
         region = gpd.GeoDataFrame()
         if "region" in self.geoms:
             region = self.geoms["region"]
-        elif "msk" in self.grid and np.any(self.grid["msk"] > 0):
-            da = xr.where(self.mask > 0, 1, 0).astype(np.int16)
-            da.raster.set_nodata(0)
-            region = da.raster.vectorize().dissolve()
-        elif self.reggrid is not None:
-            region = self.reggrid.empty_mask.raster.box
+        elif self.grid_type == "regular":
+            if "msk" in self.grid and np.any(self.grid["msk"] > 0):
+                da = xr.where(self.mask > 0, 1, 0).astype(np.int16)
+                da.raster.set_nodata(0)
+                region = da.raster.vectorize().dissolve()
+            elif self.reggrid is not None:
+                region = self.reggrid.empty_mask.raster.box
+        elif self.grid_type == "quadtree":
+            region = self.quadtree.exterior
         return region
 
+    @property
+    def bounds(self) -> List[float]:
+        """Returns the bounding box of the model grid."""
+        if self.grid_type == "regular":
+            return self.mask.raster.bounds
+        elif self.grid_type == "quadtree":
+            return self.mask.ugrid.total_bounds
+
+    @property
+    def bbox(self) -> tuple:
+        """Returns the bounding box in WGS 84 of the model grid."""
+        if self.grid_type == "regular":
+            return self.mask.raster.transform_bounds(4326)
+        elif self.grid_type == "quadtree":
+            return self.mask.ugrid.to_crs(4326).ugrid.total_bounds
+
     @property
     def crs(self) -> CRS | None:
         """Returns the model crs"""
         if self.grid_type == "regular":
             return self.reggrid.crs
         elif self.grid_type == "quadtree":
-            return self.quadtree.crs
+            return self.quadtree.data.grid.crs
 
     def set_crs(self, crs: Any) -> None:
         """Sets the model crs"""
@@ -177,7 +214,7 @@ def set_crs(self, crs: Any) -> None:
             self.reggrid.crs = CRS.from_user_input(crs)
             self.grid.raster.set_crs(self.reggrid.crs)
         elif self.grid_type == "quadtree":
-            self.quadtree.crs = CRS.from_user_input(crs)
+            self.quadtree.data.grid.set_crs(CRS.from_user_input(crs))
 
     def setup_grid(
         self,
@@ -890,7 +927,7 @@ def setup_river_inflow(
         if hydrography is not None:
             ds = self.data_catalog.get_rasterdataset(
                 hydrography,
-                bbox=self.mask.raster.transform_bounds(4326),
+                bbox=self.bbox,
                 variables=["uparea", "flwdir"],
                 buffer=5,
             )
@@ -1027,7 +1064,7 @@ def setup_river_outflow(
         if hydrography is not None:
             ds = self.data_catalog.get_rasterdataset(
                 hydrography,
-                bbox=self.mask.raster.transform_bounds(4326),
+                bbox=self.bbox,
                 variables=["uparea", "flwdir"],
                 buffer=5,
             )
@@ -1136,7 +1173,7 @@ def setup_constant_infiltration(
         if qinf is not None:
             da_inf = self.data_catalog.get_rasterdataset(
                 qinf,
-                bbox=self.mask.raster.transform_bounds(4326),
+                bbox=self.bbox,
                 buffer=10,
             )
         elif lulc is not None:
@@ -1147,7 +1184,7 @@ def setup_constant_infiltration(
                 )
             da_lulc = self.data_catalog.get_rasterdataset(
                 lulc,
-                bbox=self.mask.raster.transform_bounds(4326),
+                bbox=self.bbox,
                 buffer=10,
                 variables=["lulc"],
             )
@@ -1207,9 +1244,7 @@ def setup_cn_infiltration(self, cn, antecedent_moisture="avg", reproj_method="me
             By default 'med'. For more information see, :py:meth:`hydromt.raster.RasterDataArray.reproject_like`
         """
         # get data
-        da_org = self.data_catalog.get_rasterdataset(
-            cn, bbox=self.mask.raster.transform_bounds(4326), buffer=10
-        )
+        da_org = self.data_catalog.get_rasterdataset(cn, bbox=self.bbox, buffer=10)
         # read variable
         v = "cn"
         if antecedent_moisture:
@@ -1258,14 +1293,10 @@ def setup_cn_infiltration_with_ks(
 
         # Read the datafiles
         da_landuse = self.data_catalog.get_rasterdataset(
-            lulc, bbox=self.mask.raster.transform_bounds(4326), buffer=10
-        )
-        da_HSG = self.data_catalog.get_rasterdataset(
-            hsg, bbox=self.mask.raster.transform_bounds(4326), buffer=10
-        )
-        da_Ksat = self.data_catalog.get_rasterdataset(
-            ksat, bbox=self.mask.raster.transform_bounds(4326), buffer=10
+            lulc, bbox=self.bbox, buffer=10
         )
+        da_HSG = self.data_catalog.get_rasterdataset(hsg, bbox=self.bbox, buffer=10)
+        da_Ksat = self.data_catalog.get_rasterdataset(ksat, bbox=self.bbox, buffer=10)
         df_map = self.data_catalog.get_dataframe(reclass_table, index_col=0)
 
         # Define outputs
@@ -1341,7 +1372,7 @@ def setup_cn_infiltration_with_ks(
                 da_ks[sn, sm] = da_ks_block
 
         # Done
-        self.logger.info(f"Done with determination of values (in blocks).")
+        self.logger.info("Done with determination of values (in blocks).")
 
         # Specify the effective soil retention (seff)
         da_seff = da_smax
@@ -1464,7 +1495,7 @@ def setup_observation_points(
         if merge and name in self.geoms:
             gdf0 = self._geoms.pop(name)
             gdf_obs = gpd.GeoDataFrame(pd.concat([gdf_obs, gdf0], ignore_index=True))
-            self.logger.info(f"Adding new observation points to existing ones.")
+            self.logger.info("Adding new observation points to existing ones.")
 
         self.set_geoms(gdf_obs, name)
         self.set_config(f"{name}file", f"sfincs.{name}")
@@ -1507,7 +1538,7 @@ def setup_observation_lines(
         if merge and name in self.geoms:
             gdf0 = self._geoms.pop(name)
             gdf_obs = gpd.GeoDataFrame(pd.concat([gdf_obs, gdf0], ignore_index=True))
-            self.logger.info(f"Adding new observation lines to existing ones.")
+            self.logger.info("Adding new observation lines to existing ones.")
 
         self.set_geoms(gdf_obs, name)
         self.set_config(f"{name}file", f"sfincs.{name}")
@@ -1823,7 +1854,7 @@ def set_forcing_1d(
             if not set(gdf_locs.index) == set(df_ts.columns):
                 gdf_locs = gdf_locs.set_index(df_ts.columns)
                 self.logger.info(
-                    f"No matching index column found in gdf_locs; assuming the order is correct"
+                    "No matching index column found in gdf_locs; assuming the order is correct"
                 )
         # merge with existing data
         if name in self.forcing and merge:
@@ -1914,7 +1945,7 @@ def setup_waterlevel_forcing(
         gdf_locs, df_ts = None, None
         tstart, tstop = self.get_model_time()  # model time
         # buffer around msk==2 values
-        if np.any(self.mask == 2):
+        if np.any(self.mask == 2) and not self.grid_type == "quadtree":
             region = self.mask.where(self.mask == 2, 0).raster.vectorize()
         else:
             region = self.region
@@ -1964,7 +1995,7 @@ def setup_waterlevel_forcing(
             else:
                 da_offset = self.data_catalog.get_rasterdataset(
                     offset,
-                    bbox=self.mask.raster.transform_bounds(4326),
+                    bbox=self.bbox,
                     buffer=5,
                 )
                 offset_pnts = da_offset.raster.sample(gdf_locs)
@@ -2179,7 +2210,7 @@ def setup_discharge_forcing_from_grid(
         # read data
         ds = self.data_catalog.get_rasterdataset(
             discharge,
-            bbox=self.mask.raster.transform_bounds(4326),
+            bbox=self.bbox,
             buffer=2,
             time_tuple=self.get_model_time(),  # model time
             variables=["discharge"],
@@ -2188,7 +2219,7 @@ def setup_discharge_forcing_from_grid(
         if uparea is not None and "uparea" in gdf.columns:
             da_upa = self.data_catalog.get_rasterdataset(
                 uparea,
-                bbox=self.mask.raster.transform_bounds(4326),
+                bbox=self.bbox,
                 buffer=2,
                 variables=["uparea"],
             )
@@ -2250,7 +2281,7 @@ def setup_precip_forcing_from_grid(
         # get data for model domain and config time range
         precip = self.data_catalog.get_rasterdataset(
             precip,
-            bbox=self.mask.raster.transform_bounds(4326),
+            bbox=self.bbox,
             buffer=2,
             time_tuple=self.get_model_time(),
             variables=["precip"],
@@ -2745,11 +2776,17 @@ def plot_basemap(
         if isinstance(variable, xr.DataArray):
             ds = variable.to_dataset()
             variable = variable.name
+        elif isinstance(variable, xu.UgridDataArray):
+            ds = variable.to_dataset()
+            variable = variable.name
         elif variable.startswith("subgrid.") and self.subgrid is not None:
             ds = self.subgrid.copy()
             variable = variable.replace("subgrid.", "")
         else:
-            ds = self.grid.copy()
+            if self.grid_type == "regular":
+                ds = self.grid.copy()
+            elif self.grid_type == "quadtree":
+                ds = self.quadtree.data.copy()
             if "msk" not in ds:
                 ds["msk"] = self.mask
 
@@ -2767,6 +2804,7 @@ def plot_basemap(
             geom_names=geom_names,
             geom_kwargs=geom_kwargs,
             legend_kwargs=legend_kwargs,
+            logger=self.logger,
             **kwargs,
         )
 
@@ -2785,7 +2823,16 @@ def read(self, epsg: int = None):
         self.read_config(epsg=epsg)
         if epsg is None and "epsg" not in self.config:
             raise ValueError("Please specify epsg to read this model")
-        self.read_grid()
+        if self.grid_type == "regular":
+            self.read_grid()
+        elif self.grid_type == "quadtree":
+            fn = self.get_config("qtrfile", fallback="sfincs.nc", abs_path=True)
+            if not isfile(fn):
+                raise IOError(f".nc path {fn} does not exist")
+            self.quadtree.read(file_name=fn)
+            # remove grid from api and model
+            self._API.pop("grid", None)
+            self._grid = None
         self.read_subgrid()
         self.read_geoms()
         self.read_forcing()
@@ -2795,7 +2842,11 @@ def write(self):
         """Write the complete model schematization and configuration to file."""
         self.logger.info(f"Writing model data to {self.root}")
         # TODO - add check for subgrid & quadtree > give flags to self.write_grid() and self.write_config()
-        self.write_grid()
+        if self.grid_type == "regular":
+            self.write_grid()
+        elif self.grid_type == "quadtree":
+            fn = self.get_config("qtrfile", abs_path=True)
+            self.quadtree.write(file_name=fn)
         self.write_subgrid()
         self.write_geoms()
         self.write_forcing()
@@ -2851,6 +2902,17 @@ def read_grid(self, data_vars: Union[List, str] = None) -> None:
                 ds.raster.set_crs(epsg)
             self.set_grid(ds)
 
+            # TODO - fix this properly; but to create overlays in GUIs,
+            # we always convert regular grids to a UgridDataArray
+            self.quadtree = QuadtreeGrid(logger=self.logger)
+            if self.config.get("rotation", 0) != 0:  # This is a rotated regular grid
+                self.quadtree.data = UgridDataArray.from_structured(
+                    self.mask, "xc", "yc"
+                )
+            else:
+                self.quadtree.data = UgridDataArray.from_structured(self.mask)
+            self.quadtree.data.grid.set_crs(self.crs)
+
             # keep some metadata maps from gis directory
             fns = glob.glob(join(self.root, "gis", "*.tif"))
             fns = [
@@ -2922,29 +2984,30 @@ def read_subgrid(self):
                 self.logger.warning(f"sbgfile not found at {fn}")
                 return
 
-            # re-initialize subgrid (different variables for old/new version)
-            # TODO: come up with a better way to handle this
-            self.reggrid.subgrid = SubgridTableRegular()
-            self.subgrid = xr.Dataset()
-
-            # read subgrid file
-            if fn.parts[-1].endswith(".sbg"):  # read binary file
-                self.reggrid.subgrid.read_binary(file_name=fn, mask=self.mask)
-            else:  # read netcdf file
-                self.reggrid.subgrid.read(file_name=fn, mask=self.mask)
-            self.subgrid = self.reggrid.subgrid.to_xarray(
-                dims=self.mask.raster.dims, coords=self.mask.raster.coords
-            )
+            if self.grid_type == "regular":
+                # TODO: come up with a better way to handle this
+                self.reggrid.subgrid = SubgridTableRegular()
+                self.subgrid = xr.Dataset()
+
+                # read subgrid file
+                if fn.parts[-1].endswith(".sbg"):  # read binary file
+                    self.reggrid.subgrid.read_binary(file_name=fn, mask=self.mask)
+                else:  # read netcdf file
+                    self.reggrid.subgrid.read(file_name=fn, mask=self.mask)
+                self.subgrid = self.reggrid.subgrid.to_xarray(
+                    dims=self.mask.raster.dims, coords=self.mask.raster.coords
+                )
+            else:
+                self.quadtree.subgrid.read(file_name=fn)
 
     def write_subgrid(self):
         """Write SFINCS subgrid file."""
         self._assert_write_mode
 
-        if self.subgrid:
+        if self.grid_type == "regular" and self.subgrid:
             if "sbgfile" not in self.config:
                 # apparently no subgrid was read, so set default filename
                 self.set_config("sbgfile", "sfincs_subgrid.nc")
-
             fn = self.get_config("sbgfile", abs_path=True)
             if fn.parts[-1].endswith(".sbg"):
                 # write binary file
@@ -2952,6 +3015,13 @@ def write_subgrid(self):
             else:
                 # write netcdf file
                 self.reggrid.subgrid.write(file_name=fn, mask=self.mask)
+        elif self.grid_type == "quadtree" and self.quadtree.subgrid.data is not None:
+            if "sbgfile" not in self.config:
+                # apparently no subgrid was read, so set default filename
+                self.set_config("sbgfile", "sfincs_subgrid.nc")
+
+            fn = self.get_config("sbgfile", abs_path=True)
+            self.quadtree.subgrid.write(file_name=fn)
 
     def read_geoms(self):
         """Read geometry files and save to `geoms` attribute.
@@ -3338,16 +3408,41 @@ def read_results(
         if not isabs(fn_map):
             fn_map = join(self.root, fn_map)
         if isfile(fn_map):
-            ds_face, ds_edge = utils.read_sfincs_map_results(
-                fn_map,
-                ds_like=self.grid,  # TODO: fix for quadtree
-                drop=drop,
-                logger=self.logger,
-                **kwargs,
-            )
-            # save as dict of DataArray
-            self.set_results(ds_face, split_dataset=True)
-            self.set_results(ds_edge, split_dataset=True)
+            if self.grid_type is None:
+                self.logger.warning(
+                    "Grid type not known, trying to read from config file. "
+                )
+                self.read_config()
+            if self.grid_type == "regular":
+                ds_face, ds_edge = utils.read_sfincs_map_results(
+                    fn_map,
+                    ds_like=self.grid,  # TODO: fix for quadtree
+                    drop=drop,
+                    logger=self.logger,
+                    **kwargs,
+                )
+                # save as dict of DataArray
+                self.set_results(ds_face, split_dataset=True)
+                self.set_results(ds_edge, split_dataset=True)
+            elif self.grid_type == "quadtree":
+                dsu = xu.open_dataset(
+                    fn_map,
+                    chunks={"time": chunksize},
+                )
+
+                # set coords
+                dsu = dsu.set_coords(["mesh2d_node_x", "mesh2d_node_y"])
+                # get crs variable, drop it and set it correctly
+                crs = dsu["crs"].values
+                dsu.drop_vars("crs")
+                dsu.grid.set_crs(CRS.from_user_input(crs))
+
+                # NOTE the set_results of the model api doesnt support ugrid
+                self._initialize_results()
+                for name in dsu.variables:
+                    if name in self._results:
+                        self.logger.warning(f"Replacing result: {name}")
+                    self._results[name] = dsu[name]
 
         if not isabs(fn_his):
             fn_his = join(self.root, fn_his)
@@ -3568,9 +3663,8 @@ def update_grid_from_config(self):
                 rotation=self.config.get("rotation", 0),
                 epsg=self.config.get("epsg"),
             )
-        else:
-            raise not NotImplementedError("Quadtree grid not implemented yet")
-            # self.quadtree = QuadtreeGrid()
+        elif self.grid_type == "quadtree":
+            self.quadtree = QuadtreeGrid(logger=self.logger)
 
     def get_model_time(self):
         """Return (tstart, tstop) tuple with parsed model start and end time"""
@@ -3607,7 +3701,7 @@ def _parse_datasets_dep(self, datasets_dep, res):
                 try:
                     da_elv = self.data_catalog.get_rasterdataset(
                         dataset.get("elevtn", dataset.get("da")),
-                        bbox=self.mask.raster.transform_bounds(4326),
+                        bbox=self.bbox,
                         buffer=10,
                         variables=["elevtn"],
                         zoom_level=(res, "meter"),
@@ -3628,7 +3722,7 @@ def _parse_datasets_dep(self, datasets_dep, res):
             if "offset" in dataset and not isinstance(dataset["offset"], (float, int)):
                 da_offset = self.data_catalog.get_rasterdataset(
                     dataset.get("offset"),
-                    bbox=self.mask.raster.transform_bounds(4326),
+                    bbox=self.bbox,
                     buffer=10,
                 )
                 dd.update({"offset": da_offset})
@@ -3637,7 +3731,7 @@ def _parse_datasets_dep(self, datasets_dep, res):
             if "mask" in dataset:
                 gdf_valid = self.data_catalog.get_geodataframe(
                     dataset.get("mask"),
-                    bbox=self.mask.raster.transform_bounds(4326),
+                    bbox=self.bbox,
                 )
                 dd.update({"gdf_valid": gdf_valid})
 
@@ -3679,7 +3773,7 @@ def _parse_datasets_rgh(self, datasets_rgh):
             if "manning" in dataset or "da" in dataset:
                 da_man = self.data_catalog.get_rasterdataset(
                     dataset.get("manning", dataset.get("da")),
-                    bbox=self.mask.raster.transform_bounds(4326),
+                    bbox=self.bbox,
                     buffer=10,
                 )
                 dd.update({"da": da_man})
@@ -3691,11 +3785,11 @@ def _parse_datasets_rgh(self, datasets_rgh):
                     reclass_table = join(DATADIR, "lulc", f"{lulc}_mapping.csv")
                 if reclass_table is None:
                     raise IOError(
-                        f"Manning roughness 'reclass_table' csv file must be provided"
+                        "Manning roughness 'reclass_table' csv file must be provided"
                     )
                 da_lulc = self.data_catalog.get_rasterdataset(
                     lulc,
-                    bbox=self.mask.raster.transform_bounds(4326),
+                    bbox=self.bbox,
                     buffer=10,
                     variables=["lulc"],
                 )
@@ -3710,7 +3804,7 @@ def _parse_datasets_rgh(self, datasets_rgh):
             if "mask" in dataset:
                 gdf_valid = self.data_catalog.get_geodataframe(
                     dataset.get("mask"),
-                    bbox=self.mask.raster.transform_bounds(4326),
+                    bbox=self.bbox,
                 )
                 dd.update({"gdf_valid": gdf_valid})
 
@@ -3759,7 +3853,7 @@ def _parse_datasets_riv(self, datasets_riv):
                 else:
                     gdf_riv = self.data_catalog.get_geodataframe(
                         rivers,
-                        geom=self.mask.raster.box,
+                        bbox=self.bbox,
                         buffer=1e3,  # 1km
                     ).to_crs(self.crs)
                 # update missing attributes based on global values
@@ -3776,7 +3870,7 @@ def _parse_datasets_riv(self, datasets_riv):
             if "point_zb" in dataset:
                 gdf_zb = self.data_catalog.get_geodataframe(
                     dataset.get("point_zb"),
-                    geom=self.mask.raster.box,
+                    bbox=self.bbox,
                 )
                 dd.update({"gdf_zb": gdf_zb})
 
@@ -3793,7 +3887,7 @@ def _parse_datasets_riv(self, datasets_riv):
             if "mask" in dataset:
                 gdf_riv_mask = self.data_catalog.get_geodataframe(
                     dataset.get("mask"),
-                    geom=self.mask.raster.box,
+                    bbox=self.bbox,
                 )
                 dd.update({"gdf_riv_mask": gdf_riv_mask})
             elif "rivwth" not in gdf_riv:
diff --git a/hydromt_sfincs/subgrid.py b/hydromt_sfincs/subgrid.py
index b1e0a2ad..baf1bc7d 100644
--- a/hydromt_sfincs/subgrid.py
+++ b/hydromt_sfincs/subgrid.py
@@ -816,6 +816,40 @@ def from_xarray(self, ds_sbg):
             setattr(self, name, ds_sbg[name].values)
 
 
+class SubgridTableQuadtree:
+    # This code is still slow as it does not use numba
+
+    def __init__(self, version=0):
+        # A quadtree subgrid table contains data for EACH cell, u and v point in the quadtree mesh,
+        # regardless of the mask value!
+        self.version = version
+        self.data = None
+
+    def read(self, file_name):
+        """Read XArray dataset from netcdf file"""
+
+        if not os.path.isfile(file_name):
+            logger.info("File " + file_name + " does not exist!")
+            return
+
+        # Read from netcdf file with xarray
+        ds = xr.open_dataset(file_name)
+        # Transpose to ensure bins is first dimension (convert from FORTRAN convention in SFINCS to Python)
+        ds = ds.transpose("levels", "npuv", "np")
+        ds.close()  # Should this be closed ?
+        self.data = ds
+
+    def write(self, file_name):
+        """Write XArray dataset to netcdf file"""
+        # ensure levels is last dimension to match the FORTRAN convention in SFINCS
+        ds = self.data.transpose("npuv", "np", "levels")
+
+        # fix names to match SFINCS convention
+        # ds = ds.rename_vars({"uv_navg": "uv_navg_w", "uv_ffit": "uv_fnfit"})
+
+        ds.to_netcdf(file_name)
+
+
 @njit
 def process_tile_regular(
     mask,
@@ -1082,7 +1116,7 @@ def subgrid_q_table(
     zmax_b = np.max(dd_b)  # Maximum elevation side B
 
     zmin = max(zmin_a, zmin_b) + huthresh  # Minimum elevation of uv point
-    zmax = max(zmax_a, zmax_b)  # Maximum elevation of uv point
+    zmax = max(zmax_a, zmax_b) + huthresh  # Maximum elevation of uv point
 
     # Make sure zmax is always a bit higher than zmin
     if zmax < zmin + 0.001:
@@ -1102,48 +1136,55 @@ def subgrid_q_table(
 
         h = np.maximum(zbin - elevation, 0.0)  # water depth in each pixel
 
-        pwet[ibin] = (zbin - elevation > -1.0e-6).sum() / n
-
         # Side A
-        h_a = np.maximum(
-            zbin - dd_a, 0.0
-        )  # Depth of all pixels (but set min pixel height to zbot). Can be negative, but not zero (because zmin = zbot + huthresh, so there must be pixels below zb).
+        # Depth of all pixels (but set min pixel height to zbot). Can be negative, but not zero (because zmin = zbot + huthresh, so there must be pixels below zb).
+        h_a = np.maximum(zbin - dd_a, 0.0)
         q_a = h_a ** (5.0 / 3.0) / manning_a  # Determine 'flux' for each pixel
         q_a = np.mean(q_a)  # Grid-average flux through all the pixels
         h_a = np.mean(h_a)  # Grid-average depth through all the pixels
 
         # Side B
-        h_b = np.maximum(
-            zbin - dd_b, 0.0
-        )  # Depth of all pixels (but set min pixel height to zbot). Can be negative, but not zero (because zmin = zbot + huthresh, so there must be pixels below zb).
+        # Depth of all pixels (but set min pixel height to zbot). Can be negative, but not zero (because zmin = zbot + huthresh, so there must be pixels below zb).
+        h_b = np.maximum(zbin - dd_b, 0.0)
         q_b = h_b ** (5.0 / 3.0) / manning_b  # Determine 'flux' for each pixel
         q_b = np.mean(q_b)  # Grid-average flux through all the pixels
         h_b = np.mean(h_b)  # Grid-average depth through all the pixels
 
         # Compute q and h
-        q_all = np.mean(
-            h ** (5.0 / 3.0) / manning
-        )  # Determine grid average 'flux' for each pixel
+        # Determine grid average 'flux' for each pixel
+        q_all = np.mean(h ** (5.0 / 3.0) / manning)
         h_all = np.mean(h)  # grid averaged depth of A and B combined
         q_min = np.minimum(q_a, q_b)
         h_min = np.minimum(h_a, h_b)
 
         if option == 1:
             # Use old 1 option (weighted average of q_ab and q_all) option (min at bottom bin, mean at top bin)
-            w = (ibin) / (
-                nlevels - 1
-            )  # Weight (increase from 0 to 1 from bottom to top bin)
+            # Weight (increase from 0 to 1 from bottom to top bin)
+            w = (ibin) / (nlevels - 1)
             q = (1.0 - w) * q_min + w * q_all  # Weighted average of q_min and q_all
             hmean = h_all
 
+            # Wet fraction
+            pwet[ibin] = (zbin > elevation + huthresh).sum() / n
+
         elif option == 2:
-            # Use newer 2 option (minimum of q_a an q_b, minimum of h_a and h_b increasing to h_all, using pwet for weighting) option
-            pwet_a = (zbin - dd_a > -1.0e-6).sum() / (n / 2)
-            pwet_b = (zbin - dd_b > -1.0e-6).sum() / (n / 2)
+            # We want to make sure that, in the first layer, hmean does not exceed huthresh.
+            # This is done by making sure that the wet fraction is 0.0 in the first level on the shallowest side (i.e. if ibin==0, pwet_a or pwet_b must be 0.0).
+            # As a result, the weight w will be 0.0 in the first level on the shallowest side.
+            # At the bottom level (i.e. ibin is 0), the grid-averaged depth h_min is typically huthresh / (n / 2), assuming there is one unique pixel that is the lowest.
+            if ibin == 0:
+                # Ensure that either pwet_a or pwet_b is 0.0
+                pwet_a = (zbin > dd_a + huthresh + 1.0e-4).sum() / (n / 2)
+                pwet_b = (zbin > dd_b + huthresh + 1.0e-4).sum() / (n / 2)
+            else:
+                # Ensure that both pwet_a and pwet_b are 1.0 at the top level, so that the weight w is 1.0 and pwet[ibin] is 1.0
+                pwet_a = (zbin > dd_a + huthresh - 1.0e-4).sum() / (n / 2)
+                pwet_b = (zbin > dd_b + huthresh - 1.0e-4).sum() / (n / 2)
             # Weight increases linearly from 0 to 1 from bottom to top bin use percentage wet in sides A and B
-            w = 2 * np.minimum(pwet_a, pwet_b) / (pwet_a + pwet_b)
+            w = 2 * np.minimum(pwet_a, pwet_b) / max(pwet_a + pwet_b, 1.0e-9)
             q = (1.0 - w) * q_min + w * q_all  # Weighted average of q_min and q_all
             hmean = (1.0 - w) * h_min + w * h_all  # Weighted average of h_min and h_all
+            pwet[ibin] = 0.5 * (pwet_a + pwet_b)  # Combined pwet_a and pwet_b
 
         havg[ibin] = hmean  # conveyance depth
         nrep[ibin] = hmean ** (5.0 / 3.0) / q  # Representative n for qmean and hmean
@@ -1155,9 +1196,8 @@ def subgrid_q_table(
 
     # Determine nfit at zfit
     zfit = zmax + zmax - zmin
-    hfit = (
-        havg_top + zmax - zmin
-    )  # mean water depth in cell as computed in SFINCS (assuming linear relation between water level and water depth above zmax)
+    # mean water depth in cell as computed in SFINCS (assuming linear relation between water level and water depth above zmax)
+    hfit = havg_top + zmax - zmin
 
     # Compute q and navg
     h = np.maximum(zfit - elevation, 0.0)  # water depth in each pixel
diff --git a/hydromt_sfincs/utils.py b/hydromt_sfincs/utils.py
index 84fb9d1e..1d85451d 100644
--- a/hydromt_sfincs/utils.py
+++ b/hydromt_sfincs/utils.py
@@ -19,6 +19,7 @@
 from rasterio.rio.overview import get_maximum_overview_level
 from rasterio.windows import Window
 import xarray as xr
+import xugrid as xu
 from hydromt.io import write_xy
 from pyproj.crs.crs import CRS
 from shapely.geometry import LineString, Polygon
@@ -853,7 +854,7 @@ def read_sfincs_his_results(
 
 
 def downscale_floodmap(
-    zsmax: xr.DataArray,
+    zsmax: Union[xr.DataArray, xu.UgridDataArray],
     dep: Union[Path, str, xr.DataArray],
     hmin: float = 0.05,
     gdf_mask: gpd.GeoDataFrame = None,
@@ -900,8 +901,13 @@ def downscale_floodmap(
     hydromt.raster.RasterDataArray.to_raster
     """
     # get maximum water level
-    timedim = set(zsmax.dims) - set(zsmax.raster.dims)
+    if isinstance(zsmax, xu.UgridDataArray):
+        timedim = set(zsmax.dims) - set(zsmax.ugrid.grid.dims)
+    else:
+        timedim = set(zsmax.dims) - set(zsmax.raster.dims)
     if timedim:
+        logger.info(f"Multiple values present in {timedim} dimension.")
+        logger.info(f"Downscaling floodmap for maximum water level over {timedim}.")
         zsmax = zsmax.max(timedim)
 
     # Hydromt expects a string so if a Path is provided, convert to str
@@ -1005,21 +1011,44 @@ def downscale_floodmap(
                     if np.all(np.isnan(block_data)):
                         continue
 
-                    # TODO directly use the rasterio warp method rather than the raster.reproject see PR #145
-                    # Convert row and column indices to pixel coordinates
-                    cols, rows = np.indices((bm1 - bm0, bn1 - bn0))
-                    x_coords, y_coords = src.transform * (cols + bm0, rows + bn0)
-
-                    # Create xarray DataArray with coordinates
-                    block_dep = xr.DataArray(
-                        block_data.squeeze().transpose(),
-                        dims=("y", "x"),
-                        coords={
-                            "yc": (("y", "x"), y_coords),
-                            "xc": (("y", "x"), x_coords),
-                        },
-                    )
-                    block_dep.raster.set_crs(src.crs)
+                    # Determine if rotation is zero
+                    if src.transform[1] == 0 and src.transform[3] == 0:  # No rotation
+                        # Compute the 1D coordinates for x and y using the affine transformation
+                        x_coords = (
+                            src.transform[2]
+                            + (np.arange(bm0, bm1) + 0.5) * src.transform[0]
+                        )
+                        y_coords = (
+                            src.transform[5]
+                            + (np.arange(bn0, bn1) + 0.5) * src.transform[4]
+                        )
+
+                        # Create xarray DataArray with coordinates
+                        block_dep = xr.DataArray(
+                            block_data.squeeze(),
+                            dims=("y", "x"),
+                            coords={
+                                "y": ("y", y_coords),
+                                "x": ("x", x_coords),
+                            },
+                        )
+                    else:
+                        # Convert row and column indices to pixel coordinates
+                        cols, rows = np.meshgrid(
+                            np.arange(bm0, bm1), np.arange(bn0, bn1)
+                        )
+                        x_coords, y_coords = src.transform * (cols + 0.5, rows + 0.5)
+
+                        # Create xarray DataArray with coordinates
+                        block_dep = xr.DataArray(
+                            block_data.squeeze(),
+                            dims=("y", "x"),
+                            coords={
+                                "yc": (("y", "x"), y_coords),
+                                "xc": (("y", "x"), x_coords),
+                            },
+                        )
+                    block_dep.raster.set_crs(src.crs.to_epsg())
 
                     block_hmax = _downscale_floodmap_da(
                         zsmax=zsmax,
@@ -1031,7 +1060,7 @@ def downscale_floodmap(
 
                     with rasterio.open(floodmap_fn, "r+") as fm_tif:
                         fm_tif.write(
-                            np.transpose(block_hmax.values),
+                            block_hmax.values,
                             window=window,
                             indexes=1,
                         )
@@ -1149,7 +1178,7 @@ def build_overviews(
 
 
 def _downscale_floodmap_da(
-    zsmax: xr.DataArray,
+    zsmax: Union[xr.DataArray, xu.UgridDataArray],
     dep: xr.DataArray,
     hmin: float = 0.05,
     gdf_mask: gpd.GeoDataFrame = None,
@@ -1171,7 +1200,22 @@ def _downscale_floodmap_da(
     """
 
     # interpolate zsmax to dep grid
-    zsmax = zsmax.raster.reproject_like(dep, method=reproj_method)
+    if isinstance(zsmax, xr.DataArray):
+        zsmax = zsmax.raster.reproject_like(dep, method=reproj_method)
+
+    elif isinstance(zsmax, xu.UgridDataArray):
+        # if non-rotated grid, use xugrid rasterize_like
+        if dep.raster.transform[1] == 0 and dep.raster.transform[3] == 0:
+            zsmax = zsmax.ugrid.rasterize_like(dep)
+        # if rotated grid, use xugrid regridder
+        else:
+            # need to convert dep to unstructured to enable xugrid regridder
+            uda_dep = xu.UgridDataArray.from_structured(dep, "xc", "yc")
+            regridder = xu.CentroidLocatorRegridder(source=zsmax, target=uda_dep)
+            result = regridder.regrid(zsmax)
+            # map back to structured
+            zsmax = dep.copy(data=result.values.reshape(dep.shape))
+
     zsmax = zsmax.raster.mask_nodata()  # make sure nodata is nan
 
     # get flood depth
diff --git a/pyproject.toml b/pyproject.toml
index 707e6888..3a8e637c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -18,12 +18,13 @@ dependencies = [
     "numpy<2.0", # temp pin until bottleneck release v1.4
     "pandas",
     "pillow",
-    "pyflwdir>=0.5.5",
+    "pyflwdir>=0.5.5, <0.5.9",
     "pyproj",
     "rasterio",
     "scipy",
     "shapely",
     "xarray",
+    "xugrid>=0.12, <1.0",
 ]
 requires-python = ">=3.9" # fix tests to support older versions
 readme = "README.rst"
@@ -59,6 +60,9 @@ examples = [
     "matplotlib",
     "pillow",
 ]
+quadtree = [
+    "datashader",
+]
 
 full = ["hydromt_sfincs[test, doc, examples]"]
 
diff --git a/tests/data/sfincs_test_quadtree/gis/bnd.geojson b/tests/data/sfincs_test_quadtree/gis/bnd.geojson
new file mode 100644
index 00000000..e0479012
--- /dev/null
+++ b/tests/data/sfincs_test_quadtree/gis/bnd.geojson
@@ -0,0 +1,9 @@
+{
+"type": "FeatureCollection",
+"name": "bnd",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32633" } },
+"features": [
+{ "type": "Feature", "properties": { "stations": 1 }, "geometry": { "type": "Point", "coordinates": [ 319526.0, 5041108.0 ] } },
+{ "type": "Feature", "properties": { "stations": 2 }, "geometry": { "type": "Point", "coordinates": [ 329195.0, 5046243.0 ] } }
+]
+}
diff --git a/tests/data/sfincs_test_quadtree/gis/region.geojson b/tests/data/sfincs_test_quadtree/gis/region.geojson
new file mode 100644
index 00000000..b55ccd4e
--- /dev/null
+++ b/tests/data/sfincs_test_quadtree/gis/region.geojson
@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"name": "region",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32633" } },
+"features": [
+{ "type": "Feature", "properties": { "value": 1.0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 324990.358609079034068, 5044296.783995999023318 ], [ 324967.659084092068952, 5044341.334322208538651 ], [ 324878.558431673212908, 5044295.935272234492004 ], [ 324855.858906686247792, 5044340.485598444007337 ], [ 324588.55694942973787, 5044204.288448522798717 ], [ 324565.857424442772754, 5044248.83877473231405 ], [ 324075.803836139151827, 5043999.14399987552315 ], [ 324053.104311152128503, 5044043.694326085038483 ], [ 323919.453332523873542, 5043975.595751123502851 ], [ 323942.152857510896865, 5043931.045424913987517 ], [ 323585.750247835589107, 5043749.44922501873225 ], [ 323608.449772822554223, 5043704.898898809216917 ], [ 323430.248467984842137, 5043614.100798861123621 ], [ 323407.54894299787702, 5043658.651125070638955 ], [ 323318.448290579079185, 5043613.252075096592307 ], [ 323295.748765592055861, 5043657.80240130610764 ], [ 323162.0977869638009, 5043589.70382634550333 ], [ 323139.398261976835784, 5043634.254152555018663 ], [ 322872.096304720325861, 5043498.057002632878721 ], [ 322894.795829707290977, 5043453.506676423363388 ], [ 322805.695177288434934, 5043408.10762644931674 ], [ 322782.995652301469818, 5043452.657952658832073 ], [ 322693.894999882613774, 5043407.258902684785426 ], [ 322671.195474895648658, 5043451.809228893369436 ], [ 322582.094822476850823, 5043406.410178919322789 ], [ 322559.395297489885706, 5043450.960505128838122 ], [ 322069.341709186264779, 5043201.265730272047222 ], [ 322046.642184199299663, 5043245.816056481562555 ], [ 321957.54153178044362, 5043200.41700650844723 ], [ 321934.842006793420296, 5043244.967332717962563 ], [ 321845.74135437462246, 5043199.568282743915915 ], [ 321823.041829387657344, 5043244.118608953431249 ], [ 321778.491503178200219, 5043221.419083965942264 ], [ 321801.191028165165335, 5043176.86875775642693 ], [ 321622.989723327511456, 5043086.070657808333635 ], [ 321645.68924831453478, 5043041.520331598818302 ], [ 321556.588595895678736, 5042996.121281625702977 ], [ 321579.288120882643852, 5042951.570955416187644 ], [ 321222.885511207336094, 5042769.974755520001054 ], [ 321200.185986220370978, 5042814.525081729516387 ], [ 321066.535007592116017, 5042746.426506768912077 ], [ 321089.234532579081133, 5042701.876180559396744 ], [ 320599.180944275460206, 5042452.181405702605844 ], [ 320576.48141928849509, 5042496.731731912121177 ], [ 320398.280114450841211, 5042405.933631964027882 ], [ 320375.580589463817887, 5042450.483958173543215 ], [ 320108.278632207307965, 5042314.286808251403272 ], [ 320085.579107220342848, 5042358.837134460918605 ], [ 319996.478454801486805, 5042313.438084486871958 ], [ 319973.778929814521689, 5042357.988410696387291 ], [ 319840.127951186266728, 5042289.889835735782981 ], [ 319817.428426199301612, 5042334.440161945298314 ], [ 319683.77744757104665, 5042266.341586984694004 ], [ 319706.476972558011767, 5042221.791260775178671 ], [ 319305.524036673246883, 5042017.495535891503096 ], [ 319328.223561660211999, 5041972.945209681987762 ], [ 318927.270625775447115, 5041768.64948479924351 ], [ 318949.970150762412231, 5041724.099158589728177 ], [ 318860.869498343556188, 5041678.700108616612852 ], [ 318838.169973356591072, 5041723.250434826128185 ], [ 318749.069320937793236, 5041677.851384852081537 ], [ 318726.36979595082812, 5041722.40171106159687 ], [ 318637.269143531972077, 5041677.002661087550223 ], [ 318614.569618545006961, 5041721.552987297065556 ], [ 318213.616682660242077, 5041517.257262414321303 ], [ 318236.316207647207193, 5041472.70693620480597 ], [ 317969.014250390697271, 5041336.509786282666028 ], [ 316221.150826393452007, 5044766.884904407896101 ], [ 327358.732378748012707, 5050441.766151151619852 ], [ 329106.59580274525797, 5047011.391033026389778 ], [ 328928.394497907604091, 5046920.592933079227805 ], [ 328951.094022894569207, 5046876.042606869712472 ], [ 328772.892718056915328, 5046785.244506921619177 ], [ 328795.592243043880444, 5046740.694180712103844 ], [ 328751.041916834423319, 5046717.994655724614859 ], [ 328728.342391847458202, 5046762.544981934130192 ], [ 328639.241739428660367, 5046717.145931961014867 ], [ 328616.542214441695251, 5046761.6962581705302 ], [ 328527.441562022839207, 5046716.297208196483552 ], [ 328504.742037035874091, 5046760.847534405998886 ], [ 328192.889753569965251, 5046601.950859496369958 ], [ 328215.589278556930367, 5046557.400533286854625 ], [ 327992.837647509819362, 5046443.902908352203667 ], [ 328015.537172496784478, 5046399.352582142688334 ], [ 327926.436520077928435, 5046353.953532168641686 ], [ 327949.136045064893551, 5046309.403205959126353 ], [ 327904.585718855494633, 5046286.703680972568691 ], [ 327927.285243842517957, 5046242.153354763053358 ], [ 327838.184591423661914, 5046196.75430478900671 ], [ 327860.88411641062703, 5046152.203978579491377 ], [ 327549.03183294471819, 5045993.307303670793772 ], [ 327571.731357931683306, 5045948.756977461278439 ], [ 327438.080379303428344, 5045880.658402499742806 ], [ 327460.779904290393461, 5045836.108076291158795 ], [ 326837.075337358517572, 5045518.314726473763585 ], [ 326814.375812371552456, 5045562.865052682347596 ], [ 326769.825486162153538, 5045540.165527695789933 ], [ 326792.525011149118654, 5045495.615201487205923 ], [ 326614.323706311464775, 5045404.817101539112628 ], [ 326637.023231298429891, 5045360.266775329597294 ], [ 326503.37225267017493, 5045292.168200368992984 ], [ 326526.071777657198254, 5045247.617874159477651 ], [ 326392.420799028943293, 5045179.519299197942019 ], [ 326415.120324015908409, 5045134.968972988426685 ], [ 326058.717714340542443, 5044953.372773093171418 ], [ 326081.417239327507559, 5044908.822446883656085 ], [ 325680.464303442742676, 5044704.526722000911832 ], [ 325657.764778455777559, 5044749.077048210427165 ], [ 325390.462821199267637, 5044612.879898288287222 ], [ 325413.162346186232753, 5044568.329572078771889 ], [ 325234.961041348578874, 5044477.531472130678594 ], [ 325257.66056633554399, 5044432.981145921163261 ], [ 324990.358609079034068, 5044296.783995999023318 ] ] ] } }
+]
+}
diff --git a/tests/data/sfincs_test_quadtree/sfincs.bnd b/tests/data/sfincs_test_quadtree/sfincs.bnd
new file mode 100644
index 00000000..8c2437d0
--- /dev/null
+++ b/tests/data/sfincs_test_quadtree/sfincs.bnd
@@ -0,0 +1,2 @@
+319526.0 5041108.0
+329195.0 5046243.0
diff --git a/tests/data/sfincs_test_quadtree/sfincs.bzs b/tests/data/sfincs_test_quadtree/sfincs.bzs
new file mode 100644
index 00000000..e55f3c94
--- /dev/null
+++ b/tests/data/sfincs_test_quadtree/sfincs.bzs
@@ -0,0 +1,3 @@
+    0.0    0.00    0.25
+43200.0    0.75    1.00
+86400.0    0.00    0.25
diff --git a/tests/data/sfincs_test_quadtree/sfincs.inp b/tests/data/sfincs_test_quadtree/sfincs.inp
new file mode 100644
index 00000000..d9870053
--- /dev/null
+++ b/tests/data/sfincs_test_quadtree/sfincs.inp
@@ -0,0 +1,39 @@
+epsg                 = 32633
+latitude             = 0.0
+tref                 = 20100201 000000
+tstart               = 20100201 000000
+tstop                = 20100202 000000
+tspinup              = 3600
+dtout                = 3600.0
+dthisout             = 600.0
+dtrstout             = 0.0
+dtmaxout             = 86400
+trstout              = -999.0
+dtwnd                = 1800.0
+alpha                = 0.5
+theta                = 1.0
+huthresh             = 0.01
+manning_land         = 0.04
+manning_sea          = 0.02
+rgh_lev_land         = 0.0
+zsini                = 0.0
+rhoa                 = 1.25
+rhow                 = 1024.0
+dtmax                = 60.0
+advection            = 1
+baro                 = 1
+pavbnd               = 0
+gapres               = 101200.0
+stopdepth            = 100.0
+crsgeo               = 0
+btfilter             = 60.0
+viscosity            = 1
+bndfile              = sfincs.bnd
+bzsfile              = sfincs.bzs
+sbgfile              = sfincs_subgrid.nc
+qtrfile              = sfincs.nc
+inputformat          = bin
+outputformat         = net
+cdnrb                = 3
+cdwnd                = 0.0 28.0 50.0
+cdval                = 0.001 0.0025 0.0015
diff --git a/tests/data/sfincs_test_quadtree/sfincs.nc b/tests/data/sfincs_test_quadtree/sfincs.nc
new file mode 100644
index 00000000..1eb841f3
Binary files /dev/null and b/tests/data/sfincs_test_quadtree/sfincs.nc differ
diff --git a/tests/data/sfincs_test_quadtree/sfincs_map.nc b/tests/data/sfincs_test_quadtree/sfincs_map.nc
new file mode 100644
index 00000000..ef0cd74e
Binary files /dev/null and b/tests/data/sfincs_test_quadtree/sfincs_map.nc differ
diff --git a/tests/data/sfincs_test_quadtree/sfincs_subgrid.nc b/tests/data/sfincs_test_quadtree/sfincs_subgrid.nc
new file mode 100644
index 00000000..39388f0c
Binary files /dev/null and b/tests/data/sfincs_test_quadtree/sfincs_subgrid.nc differ
diff --git a/tests/test_1model_class.py b/tests/test_1model_class.py
index ab429604..8a9403c4 100644
--- a/tests/test_1model_class.py
+++ b/tests/test_1model_class.py
@@ -21,6 +21,9 @@
         "ini": "sfincs_test.yml",
         "example": "sfincs_test",
     },
+    "test2": {
+        "example": "sfincs_test_quadtree",
+    },
 }
 
 
@@ -32,8 +35,11 @@ def test_model_class(case):
     mod.read()
     # run test_model_api() method
     non_compliant_list = mod._test_model_api()
+    # drop non-compliant variables with "results" and "mesh" in name
+    non_compliant_list = [
+        v for v in non_compliant_list if "results" not in v and "mesh" not in v
+    ]
     assert len(non_compliant_list) == 0
-    # pass
 
 
 def test_states(mod):
@@ -395,20 +401,34 @@ def test_forcing_io(tmpdir):
     )
 
 
-def test_read_results():
-    root = TESTMODELDIR
+@pytest.mark.parametrize("case", list(_cases.keys()))
+def test_read_results(case):
+    root = join(TESTDATADIR, _cases[case]["example"])
     mod = SfincsModel(root=root, mode="r")
+    mod.read_results()
     assert all([v in mod.results for v in ["zs", "zsmax", "inp"]])
 
 
-def test_plots(mod):
-    mod.plot_forcing(fn_out="forcing.png")
-    assert isfile(join(mod.root, "figs", "forcing.png"))
-    mod.plot_basemap(fn_out="basemap.png")
-    assert isfile(join(mod.root, "figs", "basemap.png"))
-
-
 @pytest.mark.parametrize("case", list(_cases.keys()))
+def test_plots(case, tmpdir):
+    root = join(TESTDATADIR, _cases[case]["example"])
+    mod = SfincsModel(root=root, mode="r")
+    mod.read()
+    mod.plot_forcing(fn_out=join(tmpdir, "forcing.png"))
+    assert isfile(join(tmpdir, "forcing.png"))
+    fn_out = join(tmpdir, "basemap.png")
+    if case == "test2":
+        mod.plot_basemap(
+            fn_out=fn_out,
+            bmap="sat",
+            plot_bounds=False,  # does not work yet for quadtree
+        )
+    else:
+        mod.plot_basemap(fn_out=fn_out, bmap="sat")
+    assert isfile(fn_out)
+
+
+@pytest.mark.parametrize("case", list(_cases.keys())[:1])
 def test_model_build(tmpdir, case):
     # compare results with model from examples folder
     root = str(tmpdir.join(case))
diff --git a/tests/test_quadtree.py b/tests/test_quadtree.py
new file mode 100644
index 00000000..711e4552
--- /dev/null
+++ b/tests/test_quadtree.py
@@ -0,0 +1,37 @@
+from os.path import join, dirname, abspath
+import numpy as np
+from pyproj import CRS
+
+from hydromt_sfincs.quadtree import QuadtreeGrid
+
+TESTDATADIR = join(dirname(abspath(__file__)), "data")
+
+
+def test_quadtree_io(tmpdir):
+    # Initialize a QuadtreeGrid object
+    qtr = QuadtreeGrid()
+    # Read a quadtree netcdf file
+    qtr.read(join(TESTDATADIR, "sfincs_test_quadtree", "sfincs.nc"))
+    # Check the face coordinates
+    face_coordinates = qtr.face_coordinates
+    assert len(face_coordinates[0] == 4452)
+    # Check the msk variable
+    msk = qtr.data["msk"]
+    assert np.sum(msk.values) == 4298
+    # Check the crs
+    crs = qtr.crs
+    assert crs == CRS.from_epsg(32633)
+
+    # now write the quadtree to a new file
+    fn = tmpdir.join("sfincs_out.nc")
+    qtr.write(fn)
+
+    # read the new file and check the msk variable
+    qtr2 = QuadtreeGrid()
+    qtr2.read(fn)
+    # assert the crs is the same
+    assert qtr2.crs == qtr.crs
+    # assert the msk variable is the same
+    assert np.sum(qtr2.data["msk"].values) == 4298
+    # assert the dep variable is the same
+    assert np.sum(qtr.data["dep"].values) == np.sum(qtr2.data["dep"].values)