diff --git a/dfm_tools/bathymetry.py b/dfm_tools/bathymetry.py
index 371365bdb..cf7af068e 100644
--- a/dfm_tools/bathymetry.py
+++ b/dfm_tools/bathymetry.py
@@ -6,6 +6,7 @@
 """
 
 import numpy as np
+import xarray as xr
 
 
 def write_bathy_toasc(filename_asc,lon_sel_ext,lat_sel_ext,elev_sel_ext,asc_fmt='%9.2f',nodata_val=32767):
@@ -39,3 +40,71 @@ def write_bathy_toasc(filename_asc,lon_sel_ext,lat_sel_ext,elev_sel_ext,asc_fmt=
         np.savetxt(file_asc,np.flip(elev_sel_ext,axis=0),fmt=asc_fmt)
     print('...finished')
 
+
+def read_asc(file_asc:str) -> xr.Dataset:
+    """
+    Reading asc file into a xarray.Dataset
+
+    Parameters
+    ----------
+    file_asc : str
+        asc file with header ncols, nrows, xllcenter, yllcenter, cellsize, NODATA_value.
+
+    Returns
+    -------
+    ds_asc : xr.Dataset
+        xarray.Dataset with the data from the asc file as an array and 
+        the lat/lon coordinates as separate coordinate variables.
+
+    """
+    with open(file_asc) as f:
+        lines = f.readlines()
+    
+    # read header
+    header_dict = {}
+    for il, line in enumerate(lines):
+        linesplit = line.split()
+        linestart = linesplit[0]
+        linestop = linesplit[-1]
+        
+        try:
+            # try to convert it to float, this will fail for headerlines
+            # it will succeed for numeric lines, but then the loop breaks
+            float(linestart)
+            skiprows = il
+            break
+        except ValueError:
+            # convert header values to int if possible or float if not
+            try:
+                header_value = int(linestop)
+            except ValueError:
+                header_value = float(linestop)
+            header_dict[linestart] = header_value
+    
+    # read data
+    asc_np = np.loadtxt(file_asc, skiprows=skiprows, dtype=float)
+    
+    # derive x/y values and assert shape
+    num_x = header_dict['ncols']
+    num_y = header_dict['nrows']
+    start_x = header_dict['xllcenter']
+    start_y = header_dict['yllcenter']
+    step = header_dict['cellsize']
+    nodata = header_dict['NODATA_value']
+    
+    assert asc_np.shape == (num_y, num_x)
+    
+    x_vals = np.arange(0, num_x) * step + start_x
+    y_vals = np.arange(0, num_y) * step + start_y
+    
+    # flip over latitude and replace nodata with nan
+    asc_np = np.flipud(asc_np)
+    asc_np[asc_np==nodata] = np.nan
+    
+    ds_asc = xr.Dataset()
+    ds_asc['lon'] = xr.DataArray(x_vals, dims=('lon'))
+    ds_asc['lat'] = xr.DataArray(y_vals, dims=('lat'))
+    ds_asc['data'] = xr.DataArray(asc_np, dims=('lat','lon'))
+    ds_asc = ds_asc.assign_attrs(header_dict)
+    return ds_asc
+
diff --git a/docs/whats-new.md b/docs/whats-new.md
index 040077d84..8db7dcc6f 100644
--- a/docs/whats-new.md
+++ b/docs/whats-new.md
@@ -1,10 +1,11 @@
 ## UNRELEASED
 
 ### Feat
+- support for reading asc files with `dfmt.read_asc()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#571](https://github.com/Deltares/dfm_tools/pull/571)
 - added `dfmt.meshkernel_delete_withshp()` to delete parts of a mesh with a shapefile, while only reading the shapefile within the bounding box of the mesh by [@rqwang](https://github.com/rqwang) in [#548](https://github.com/Deltares/dfm_tools/pull/548) and [#566](https://github.com/Deltares/dfm_tools/pull/566)
-- improved spatial interpolation in `dfmt.interp_regularnc_to_plipoints()` by combining linear with nearest interpolation by [@veenstrajelmer] in [#561](https://github.com/Deltares/dfm_tools/pull/561)
-- added `GTSMv4.1` and `GTSMv4.1_opendap` datasets for tide interpolation with `dfmt.interpolate_tide_to_bc()` by [@veenstrajelmer] in [#544](https://github.com/Deltares/dfm_tools/pull/544)
-- support for `preprocess` argument in `dfmt.open_partitioned_dataset()`  by [@veenstrajelmer] in [#530](https://github.com/Deltares/dfm_tools/pull/530)
+- improved spatial interpolation in `dfmt.interp_regularnc_to_plipoints()` by combining linear with nearest interpolation by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#561](https://github.com/Deltares/dfm_tools/pull/561)
+- added `GTSMv4.1` and `GTSMv4.1_opendap` datasets for tide interpolation with `dfmt.interpolate_tide_to_bc()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#544](https://github.com/Deltares/dfm_tools/pull/544)
+- support for `preprocess` argument in `dfmt.open_partitioned_dataset()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#530](https://github.com/Deltares/dfm_tools/pull/530)
 
 
 ## 0.14.0 (2023-09-15)
diff --git a/tests/test_bathymetry.py b/tests/test_bathymetry.py
new file mode 100644
index 000000000..fdbaeca7b
--- /dev/null
+++ b/tests/test_bathymetry.py
@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Oct 11 10:48:06 2023
+
+@author: veenstra
+"""
+
+import os
+import pytest
+import numpy as np
+import dfm_tools as dfmt
+
+
+@pytest.mark.unittest
+def test_read_write_asc():
+    lat_range = np.arange(-8,2,0.25)
+    lon_range = np.arange(-10,1,0.25)
+    data = np.cos(lon_range) * lat_range[np.newaxis].T
+    file_asc = './dummy.asc'
+    dfmt.write_bathy_toasc(file_asc, lon_sel_ext=lon_range, lat_sel_ext=lat_range, elev_sel_ext=data, asc_fmt='%14.9f',nodata_val=-999)
+    
+    ds_asc = dfmt.read_asc(file_asc)
+    
+    assert (np.abs(ds_asc['lon'].to_numpy() - lon_range) < 1e-9).all()
+    assert (np.abs(ds_asc['lat'].to_numpy() - lat_range) < 1e-9).all()
+    assert (np.abs(ds_asc['data'].to_numpy() - data) < 1e-9).all()
+    
+    os.remove(file_asc)