Repository: trhallam/segysak (language: python)

Issue:
freeloader to use `map` instead of `replace` for performance
I am working through using `segy_freeloader` and things were going slow.
I broke down the function and noticed the `replace` function was taking time.
[This StackOverflow answer](https://stackoverflow.com/a/41987136/12127259) points out that `map` may be faster.

I ran a [timing test](https://gist.github.com/fr1ll/7673abe5b3f1b640bc6c0e0fe2feaf8c) for series of a range of sizes.
In my test, `map` is 100x faster for larger series.
![2020-10-05_timing-test-map-vs-replace](https://user-images.githubusercontent.com/29168593/95152631-24710a00-0753-11eb-9150-974f38f84db5.png)

Please let me know if my style is problematic for this PR, I am new at this.
I read the style guide and noticed that normally we should submit an issue prior to a PR,
but since this was a smallish change, thought I'd just submit the PR.




Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index ce2cc4b..7fa7ec3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,6 +24,11 @@ wheels/
 .installed.cfg
 *.egg
 MANIFEST
+
+# created docs
+docs/examples
+examples/*.ipynb
+
 # don't commit auto-generated version.py
 segysak/version.py
 
diff --git a/README.rst b/README.rst
index 810b594..bdfc9c3 100644
--- a/README.rst
+++ b/README.rst
@@ -136,4 +136,4 @@ The complete documentation for *SEGY-SAK* can be found at
 
 .. |slack| image:: https://img.shields.io/badge/slack-help-blueviolet?logo=slack
    :alt: help
-   :target: https://swung.slack.com/messages/t20-segysak/
+   :target: https://swung.slack.com/messages/segysak/
diff --git a/contributing.rst b/contributing.rst
index dfed09d..540603f 100644
--- a/contributing.rst
+++ b/contributing.rst
@@ -9,7 +9,7 @@ These are some of the many ways to contribute:
 
 * Submitting bug reports and feature requests
 * Writing tutorials or examples
-* Fixing typos and improving to the documentation
+* Fixing typos and improving the documentation
 * Writing code for everyone to use
 
 If you get stuck at any point you can create an issue on GitHub (look for the *Issues*
diff --git a/docs/examples.rst b/docs/examples.rst
index 5fab7d9..2922d48 100644
--- a/docs/examples.rst
+++ b/docs/examples.rst
@@ -9,5 +9,6 @@ Examples
     examples/example_segy_headers
     examples/example_amplitude_extraction_displays
     examples/example_extract_arbitrary_line
+    examples/example_subsurface
     examples/example_segysak_segy_vectorisation
     examples/example_segysak_dask
diff --git a/docs/requirements.txt b/docs/requirements.txt
index e689f0a..0843587 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -9,3 +9,4 @@ xarray
 numpy
 nbsphinx
 click >= 7.0, <8.0
+more-itertools
diff --git a/environment.yml b/environment.yml
index 056ef0d..41a989e 100644
--- a/environment.yml
+++ b/environment.yml
@@ -28,6 +28,7 @@ dependencies:
   - click
   - attrdict
   - pyshp
+  - more-itertools
   - pip:
       - sphinxcontrib-programoutput
       - sphinx-copybutton
diff --git a/examples/example_extract_arbitrary_line.py b/examples/example_extract_arbitrary_line.py
index d0cd514..ab6b752 100644
--- a/examples/example_extract_arbitrary_line.py
+++ b/examples/example_extract_arbitrary_line.py
@@ -24,7 +24,8 @@
 
 # %% nbsphinx="hidden"
 import warnings
-warnings.filterwarnings('ignore')
+
+warnings.filterwarnings("ignore")
 
 # %%
 from os import path
@@ -39,7 +40,13 @@
 print(f"{volve_3d_path} exists: {path.exists(volve_3d_path)}")
 
 # %%
-from segysak.segy import segy_loader, get_segy_texthead, segy_header_scan, segy_header_scrape, well_known_byte_locs
+from segysak.segy import (
+    segy_loader,
+    get_segy_texthead,
+    segy_header_scan,
+    segy_header_scrape,
+    well_known_byte_locs,
+)
 
 volve_3d = segy_loader(volve_3d_path, **well_known_byte_locs("petrel_3d"))
 
@@ -49,8 +56,14 @@
 # We define a line as lists of cdp_x & cdp_y points. These can be inside or outside of the survey, but obviously should intersect it in order to be useful.
 
 # %%
-arb_line_A = (np.array([434300, 434600, 435500, 436300]), np.array([6.4786e6, 6.4780e6, 6.4779e6, 6.4781e6]))
-arb_line_B = (np.array([434000, 434600, 435500, 436500]), np.array([6.4786e6, 6.4776e6, 6.4786e6, 6.4775e6]))
+arb_line_A = (
+    np.array([434300, 434600, 435500, 436300]),
+    np.array([6.4786e6, 6.4780e6, 6.4779e6, 6.4781e6]),
+)
+arb_line_B = (
+    np.array([434000, 434600, 435500, 436500]),
+    np.array([6.4786e6, 6.4776e6, 6.4786e6, 6.4775e6]),
+)
 
 # %% [markdown]
 # Let's see how these lines are placed relative to the survey bounds. We can see *A* is full enclosed whilst *B* has some segments outside.
@@ -62,7 +75,7 @@
 ax.legend()
 
 # %% [markdown]
-# Let's extract line A. 
+# Let's extract line A.
 #
 # We specify a `bin_spacing_hint` which is our desired bin spacing along the line. The function use this hint to calculate the closest binspacing that maintains uniform sampling.
 #
@@ -77,7 +90,7 @@
 print(f"That took {toc-tic} seconds")
 
 # %%
-plt.figure(figsize=(8,8))
+plt.figure(figsize=(8, 8))
 plt.imshow(line_A.to_array().squeeze().T, aspect="auto", cmap="RdBu", vmin=-10, vmax=10)
 
 # %% [markdown]
@@ -90,13 +103,13 @@
 print(f"That took {toc-tic} seconds")
 
 # %%
-plt.figure(figsize=(10,8))
+plt.figure(figsize=(10, 8))
 plt.imshow(line_B.to_array().squeeze().T, aspect="auto", cmap="RdBu", vmin=-10, vmax=10)
 
 # %% [markdown]
 # ## Petrel Shapefile
 #
-# We have an arbitrary line geometry defined over this small survey region stored in a shape file exported from Petrel. 
+# We have an arbitrary line geometry defined over this small survey region stored in a shape file exported from Petrel.
 #
 # Let's load that and extract an arbitrary line using segysak. We also have the seismic data extracted along that line by Petrel, so we can see how that compares.
 
@@ -108,8 +121,8 @@
 # Load the shapefile and get the list of points
 
 # %%
-sf = shapefile.Reader(path.join("data","arbitrary_line.shp"))
-line_petrel = sf.shapes() 
+sf = shapefile.Reader(path.join("data", "arbitrary_line.shp"))
+line_petrel = sf.shapes()
 print(f"shapes are type {sf.shapeType} == POLYLINEZ")
 print(f"There are {len(sf.shapes())} shapes in here")
 print(f"The line has {len(sf.shape(0).points)} points")
@@ -122,18 +135,24 @@
 
 # %%
 line_extracted_by_petrel_path = path.join("data", "volve10r12-full-twt-arb.sgy")
-print(f"{line_extracted_by_petrel_path} exists: {path.exists(line_extracted_by_petrel_path)}")
-line_extracted_by_petrel = segy_loader(line_extracted_by_petrel_path, **well_known_byte_locs("petrel_3d"))
+print(
+    f"{line_extracted_by_petrel_path} exists: {path.exists(line_extracted_by_petrel_path)}"
+)
+line_extracted_by_petrel = segy_loader(
+    line_extracted_by_petrel_path, **well_known_byte_locs("petrel_3d")
+)
 
 # %% [markdown]
 # Extract the line using segysak
 
 # %%
 tic = time()
-line_extracted_by_segysak = volve_3d.seis.interp_line(*points_from_shapefile,
-                                        bin_spacing_hint=10,
-                                        line_method='linear',
-                                        xysel_method='linear')
+line_extracted_by_segysak = volve_3d.seis.interp_line(
+    *points_from_shapefile,
+    bin_spacing_hint=10,
+    line_method="linear",
+    xysel_method="linear",
+)
 toc = time()
 print(f"That took {toc-tic} seconds")
 
@@ -141,18 +160,30 @@
 # Plot the extracted lines side by side
 
 # %%
-fig, axs = plt.subplots(1,2, figsize=(16,8))
-
-axs[0].imshow(line_extracted_by_segysak.to_array().squeeze().T, aspect="auto", cmap="RdBu", vmin=-10, vmax=10)
+fig, axs = plt.subplots(1, 2, figsize=(16, 8))
+
+axs[0].imshow(
+    line_extracted_by_segysak.to_array().squeeze().T,
+    aspect="auto",
+    cmap="RdBu",
+    vmin=-10,
+    vmax=10,
+)
 axs[0].set_title("segysak")
-axs[1].imshow(line_extracted_by_petrel.to_array().squeeze().T, aspect="auto", cmap="RdBu", vmin=-10, vmax=10)
+axs[1].imshow(
+    line_extracted_by_petrel.to_array().squeeze().T,
+    aspect="auto",
+    cmap="RdBu",
+    vmin=-10,
+    vmax=10,
+)
 axs[1].set_title("petrel")
 
 # %% [markdown]
 # Plot the geometry, trace header locations along with the volve 3d bound box, to make sure things line up.
 
 # %%
-plt.figure(figsize=(10,10))
+plt.figure(figsize=(10, 10))
 ax = volve_3d.seis.plot_bounds(ax=plt.gca())
 
 # plot path
@@ -160,9 +191,19 @@
 ax.scatter(*points_from_shapefile)
 
 # plot trace positons from extracted lines based on header
-ax.scatter(line_extracted_by_segysak.cdp_x, line_extracted_by_segysak.cdp_y, marker='x', color='k')
-ax.scatter(line_extracted_by_petrel.cdp_x, line_extracted_by_petrel.cdp_y, marker='+', color='r')
-
-ax.set_xlim(434500,435500)
+ax.scatter(
+    line_extracted_by_segysak.cdp_x,
+    line_extracted_by_segysak.cdp_y,
+    marker="x",
+    color="k",
+)
+ax.scatter(
+    line_extracted_by_petrel.cdp_x,
+    line_extracted_by_petrel.cdp_y,
+    marker="+",
+    color="r",
+)
+
+ax.set_xlim(434500, 435500)
 ax.set_ylim(6477800, 6478600)
-plt.legend(labels=['bounds','geom','corner','segysak','petrel'])
+plt.legend(labels=["bounds", "geom", "corner", "segysak", "petrel"])
diff --git a/examples/example_extract_data_on_a_horizon.py b/examples/example_extract_data_on_a_horizon.py
index 6796fdd..aeae9d9 100644
--- a/examples/example_extract_data_on_a_horizon.py
+++ b/examples/example_extract_data_on_a_horizon.py
@@ -18,16 +18,23 @@
 
 # %% nbsphinx="hidden"
 import warnings
-warnings.filterwarnings('ignore')
+
+warnings.filterwarnings("ignore")
 
 # %%
 # %load_ext autoreload
 # %autoreload 2
 
 import matplotlib.pyplot as plt
+
 # %matplotlib inline
 
-from segysak.segy import segy_loader, get_segy_texthead, segy_header_scan, segy_header_scrape
+from segysak.segy import (
+    segy_loader,
+    get_segy_texthead,
+    segy_header_scan,
+    segy_header_scrape,
+)
 
 from os import path
 
@@ -57,7 +64,7 @@
 # %%
 import pandas as pd
 
-top_hugin_df = pd.read_csv(top_hugin_path, names=["cdp_x","cdp_y","twt"], sep=' ')
+top_hugin_df = pd.read_csv(top_hugin_path, names=["cdp_x", "cdp_y", "twt"], sep=" ")
 top_hugin_df.head()
 
 # %% [markdown]
@@ -73,7 +80,9 @@
 # Alternative we can use the points to output a `xarray.Dataset` which comes with coordinates for plotting already gridded up for Pyvista.
 
 # %%
-top_hugin_ds = volve_3d.seis.surface_from_points(top_hugin_df, 'twt', right=('cdp_x', 'cdp_y'))
+top_hugin_ds = volve_3d.seis.surface_from_points(
+    top_hugin_df, "twt", right=("cdp_x", "cdp_y")
+)
 top_hugin_ds
 
 # %%
@@ -94,11 +103,15 @@
 # Extracting horizon amplitudes requires us to interpolate the cube onto the 3D horizon.
 
 # %%
-top_hugin_amp = volve_3d.data.interp({'iline':top_hugin_ds.iline, 'xline':top_hugin_ds.xline, 'twt':top_hugin_ds.twt})
+top_hugin_amp = volve_3d.data.interp(
+    {"iline": top_hugin_ds.iline, "xline": top_hugin_ds.xline, "twt": top_hugin_ds.twt}
+)
 
 # %%
-# 
+#
 fig = plt.figure(figsize=(15, 5))
-top_hugin_amp.plot(cmap='bwr')
-cs = plt.contour(top_hugin_amp.xline, top_hugin_amp.iline, top_hugin_ds.twt, levels=20, colors='grey')
-plt.clabel(cs, fontsize=14, fmt='%.1f')
+top_hugin_amp.plot(cmap="bwr")
+cs = plt.contour(
+    top_hugin_amp.xline, top_hugin_amp.iline, top_hugin_ds.twt, levels=20, colors="grey"
+)
+plt.clabel(cs, fontsize=14, fmt="%.1f")
diff --git a/examples/example_segysak_dask.py b/examples/example_segysak_dask.py
index fa0fa3a..d581807 100644
--- a/examples/example_segysak_dask.py
+++ b/examples/example_segysak_dask.py
@@ -27,13 +27,15 @@
 
 # %% nbsphinx="hidden"
 import warnings
-warnings.filterwarnings('ignore')
+
+warnings.filterwarnings("ignore")
 
 # %%
 import numpy as np
 from segysak import open_seisnc, segy
 
 import matplotlib.pyplot as plt
+
 # %matplotlib inline
 
 # %%
@@ -46,7 +48,7 @@
 # We can also scale the cluster to be a bit smaller.
 
 # %%
-client.cluster.scale(2, memory='0.5gb')
+client.cluster.scale(2, memory="0.5gb")
 client
 
 # %% [markdown]
@@ -55,9 +57,11 @@
 # If your data is in SEG-Y to use dask it must be converted to SEISNC. If you do this with the CLI it only need happen once.
 
 # %%
-segy_file = 'data/volve10r12-full-twt-sub3d.sgy'
-seisnc_file = 'data/volve10r12-full-twt-sub3d.seisnc'
-segy.segy_converter(segy_file, seisnc_file, iline=189, xline=193, cdpx=181, cdpy=185, silent=True)
+segy_file = "data/volve10r12-full-twt-sub3d.sgy"
+seisnc_file = "data/volve10r12-full-twt-sub3d.seisnc"
+segy.segy_converter(
+    segy_file, seisnc_file, iline=189, xline=193, cdpx=181, cdpy=185, silent=True
+)
 
 # %% [markdown]
 # By specifying the chunks argument to the `open_seisnc` command we can ask dask to fetch the data in chunks of size *n*. In this example the `iline` dimension will be chunked in groups of 100. The valid arguments to chunks depends on the dataset but any dimension can be used.
@@ -65,7 +69,7 @@
 # Even though the seis of the dataset is `2.14GB` it hasn't yet been loaded into memory, not will `dask` load it entirely unless the operation demands it.
 
 # %%
-seisnc = open_seisnc('data/volve10r12-full-twt-sub3d.seisnc', chunks={'iline':100})
+seisnc = open_seisnc("data/volve10r12-full-twt-sub3d.seisnc", chunks={"iline": 100})
 seisnc.seis.humanbytes
 
 # %% [markdown]
@@ -97,9 +101,9 @@
 # %%
 fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 10))
 
-iline = seisnc.sel(iline = 10100).transpose('twt', 'xline').data
-xline = seisnc.sel(xline = 2349).transpose('twt', 'iline').data
-zslice = seisnc.sel(twt = 2900, method='nearest').transpose('iline', 'xline').data
+iline = seisnc.sel(iline=10100).transpose("twt", "xline").data
+xline = seisnc.sel(xline=2349).transpose("twt", "iline").data
+zslice = seisnc.sel(twt=2900, method="nearest").transpose("iline", "xline").data
 
 q = iline.quantile([0, 0.001, 0.5, 0.999, 1]).values
 rq = np.max(np.abs([q[1], q[-2]]))
@@ -109,12 +113,12 @@
 zslice.plot(robust=True, ax=axs[0, 2])
 
 imshow_kwargs = dict(
-    cmap='seismic', aspect='auto', vmin=-rq, vmax=rq, interpolation='bicubic'
+    cmap="seismic", aspect="auto", vmin=-rq, vmax=rq, interpolation="bicubic"
 )
 
 axs[1, 0].imshow(iline.values, **imshow_kwargs)
-axs[1, 0].set_title('iline')
+axs[1, 0].set_title("iline")
 axs[1, 1].imshow(xline.values, **imshow_kwargs)
-axs[1, 1].set_title('xline')
-axs[1, 2].imshow(zslice.values, origin='lower', **imshow_kwargs)
-axs[1, 2].set_title('twt')
+axs[1, 1].set_title("xline")
+axs[1, 2].imshow(zslice.values, origin="lower", **imshow_kwargs)
+axs[1, 2].set_title("twt")
diff --git a/examples/example_subsurface.py b/examples/example_subsurface.py
new file mode 100644
index 0000000..e85b37b
--- /dev/null
+++ b/examples/example_subsurface.py
@@ -0,0 +1,45 @@
+# ---
+# jupyter:
+#   jupytext:
+#     formats: ipynb,py:percent
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.11.1
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # Using `segysak` with `subsurface`.
+
+# %%
+from segysak.segy import segy_loader
+from IPython import display
+import subsurface as ss
+import warnings
+
+warnings.filterwarnings("ignore")
+
+# %% [markdown]
+# Load in some test data.
+
+# %%
+segycube = segy_loader("data/volve10r12-full-twt-sub3d.sgy")
+
+# %%
+ss_cube = segycube.seisio.to_subsurface()
+# ss_cube is Structured Data in subsurface
+display(ss_cube)
+
+# %% [markdown]
+# Cast to a `ss.StructuredGrid` and plot using PyVista
+
+# %%
+ss_grid = ss.StructuredGrid(ss_cube)
+pvgrid = ss.visualization.to_pyvista_grid(ss_grid, "data")
+_ = ss.visualization.pv_plot([pvgrid])
+_.close()
diff --git a/segysak/_accessor.py b/segysak/_accessor.py
index 8ce155f..76ed7b0 100644
--- a/segysak/_accessor.py
+++ b/segysak/_accessor.py
@@ -48,6 +48,38 @@ def to_netcdf(self, seisnc, **kwargs):
 
         self._obj.to_netcdf(seisnc, **kwargs)
 
+    def to_subsurface(self):
+        """Convert seismic data to a subsurface StructuredData Object
+
+        Raises:
+            err: [description]
+            NotImplementedError: [description]
+
+        Returns:
+            subsurface.structs.base_structures.StructuredData: subsurface struct
+        """
+        try:
+            from subsurface.structs.base_structures import StructuredData
+        except ImportError as err:
+            print("subsurface optional dependency required")
+            raise err
+
+        if self._obj.seis.is_twt():
+            vdom = "twt"
+        else:
+            vdom = "depth"
+
+        if self._obj.seis.is_3d():
+            ds = self._obj.rename_dims({"iline": "x", "xline": "y", vdom: "z"})
+            ds["z"] = ds["z"] * -1
+            subsurface_struct = StructuredData(ds, "data")
+        elif self._obj.seis.is_2d():
+            subsurface_struct = StructuredData(self._obj, "data")
+        else:
+            raise NotImplementedError
+
+        return subsurface_struct
+
 
 def open_seisnc(seisnc, **kwargs):
     """Load from netcdf4 with seisnc specs.
@@ -318,6 +350,10 @@ def is_3d(self):
         dim_options = [
             DimensionKeyField.threed_twt,
             DimensionKeyField.threed_depth,
+            DimensionKeyField.threed_xline_twt,
+            DimensionKeyField.threed_iline_twt,
+            DimensionKeyField.threed_xline_depth,
+            DimensionKeyField.threed_iline_depth,
         ]
         invalid_dim_options = [
             DimensionKeyField.threed_ps_twt,
@@ -532,7 +568,7 @@ def calc_corner_points(self):
                 )
 
         elif self.is_2d() or self.is_2dgath():
-            cdp = DimensionKeyField.cdp_2d
+            cdp = DimensionKeyField.cdp_2d[0]
             cdps = self._obj[cdp].values
             corner_points = (cdps[0], cdps[-1])
 
@@ -620,7 +656,7 @@ def plot_bounds(self, ax=None):
         return ax
 
     def subsample_dims(self, **dim_kwargs):
-        """Return a dictionary of subsampled dims suitable for xarra.interp.
+        """Return a dictionary of subsampled dims suitable for xarray.interp.
 
         This tool halves
 
@@ -692,3 +728,37 @@ def get_affine_transform(self):
             .translate(cp_ix[0][0], cp_ix[0][1])  # move to ilxl origin
         )
         return affine_grid2loc.inverted()
+
+    def get_dead_trace_map(self, scan=None, zeros_as_nan=False):
+        """Scan the vertical axis of a volume to find traces that are all NaN
+        and return an DataArray which maps the all dead traces.
+
+        Faster scans can be performed by setting scan to an int or list of int
+        representing horizontal slice indexes to use for the scan.
+
+        Args:
+            scan (int/list of int, optional): Horizontal indexes to scan.
+                Defaults to None (scan full volume).
+            zeros_as_nan (bool, optional): Treat zeros as NaN during scan.
+
+        Returns:
+            xarray.DataArray: boolean dead trace map.
+        """
+
+        if self.is_twt():
+            vdom = CoordKeyField.twt
+        else:
+            vdom = CoordKeyField.depth
+
+        if scan:
+            nan_map = (
+                self._obj.data.isel(**{vdom: scan}).isnull().reduce(np.all, dim=vdom)
+            )
+        else:
+            nan_map = self._obj.data.isnull().reduce(np.all, dim=vdom)
+
+        if zeros_as_nan:
+            zero_map = np.abs(self._obj.data).sum(dim=vdom)
+            nan_map = xr.where(zero_map == 0.0, 1, nan_map)
+
+        return nan_map
diff --git a/segysak/_cli.py b/segysak/_cli.py
index 4b13217..4cd446d 100644
--- a/segysak/_cli.py
+++ b/segysak/_cli.py
@@ -183,7 +183,11 @@ def scrape(filename, ebcidc=False, trace_headers=False):
 @cli.command(
     help="Convert file between SEGY and NETCDF (direction is guessed or can be made explicit with the --output-type option)"
 )
-@click.argument("input-files", type=click.Path(exists=True), nargs=-1,)
+@click.argument(
+    "input-files",
+    type=click.Path(exists=True),
+    nargs=-1,
+)
 @click.option(
     "--output-file", "-o", type=click.STRING, help="Output file name", default=None
 )
@@ -220,7 +224,9 @@ def convert(
 ):
 
     if len(input_files) > 1 and output_file is not None:
-        raise ValueError("The output file option should not be used with multiple input files.")
+        raise ValueError(
+            "The output file option should not be used with multiple input files."
+        )
 
     for input_file in input_files:
         input_file = pathlib.Path(input_file)
@@ -286,7 +292,9 @@ def convert(
             click.echo(f"Converted file saved as {output_file_loc}")
             LOGGER.info(f"SEGY output written to {output_file_loc}")
         else:
-            click.echo(f"Conversion to output-type {output_type} is not implemented yet")
+            click.echo(
+                f"Conversion to output-type {output_type} is not implemented yet"
+            )
             raise SystemExit
 
 
diff --git a/segysak/_core.py b/segysak/_core.py
index 6a50451..d012f51 100644
--- a/segysak/_core.py
+++ b/segysak/_core.py
@@ -5,8 +5,7 @@
 
 
 class FrozenDict(AttrDict):
-    """A Frozen Attribute Dictionary for protected key pair mappings in SEGYSAK.
-    """
+    """A Frozen Attribute Dictionary for protected key pair mappings in SEGYSAK."""
 
     def __setitem__(self, key, value):
         raise TypeError(
diff --git a/segysak/_keyfield.py b/segysak/_keyfield.py
index 7368cab..7b77463 100644
--- a/segysak/_keyfield.py
+++ b/segysak/_keyfield.py
@@ -23,7 +23,11 @@
         cdp_3d=("iline", "xline"),
         cdp_2d=("cdp",),
         threed_twt=("iline", "xline", "twt"),
+        threed_xline_twt=("iline", "twt"),
+        threed_iline_twt=("xline", "twt"),
         threed_depth=("iline", "xline", "depth"),
+        threed_xline_depth=("iline", "depth"),
+        threed_iline_depth=("xline", "depth"),
         threed_head=("iline", "xline"),
         threed_ps_twt=("iline", "xline", "twt", "offset"),
         threed_ps_depth=("iline", "xline", "depth", "offset"),
@@ -60,6 +64,7 @@
 
 
 VerticalKeyField = FrozenDict(dict(twt="TWT", depth="DEPTH"))
+VerticalKeyDim = FrozenDict(dict(TWT="twt", DEPTH="depth"))
 
 
 VerticalUnits = FrozenDict(
diff --git a/segysak/_richstr.py b/segysak/_richstr.py
index 4af3a9b..d76b332 100644
--- a/segysak/_richstr.py
+++ b/segysak/_richstr.py
@@ -27,4 +27,3 @@ def _upgrade_txt_richstr(text):
     return type(
         "rich_texthead", (str,), dict(_repr_pretty_=_prettify, _repr_html_=_htmlify)
     )(text)
-
diff --git a/segysak/_seismic_dataset.py b/segysak/_seismic_dataset.py
index dc0bf46..0de5767 100644
--- a/segysak/_seismic_dataset.py
+++ b/segysak/_seismic_dataset.py
@@ -186,13 +186,22 @@ def create3d_dataset(
     # 1e-10 to stabilise range on rounding errors for weird floats from hypothesis
 
     vert = np.arange(
-        first_sample, first_sample + sample_rate * ns - 1e-10, sample_rate, dtype=int,
+        first_sample,
+        first_sample + sample_rate * ns - 1e-10,
+        sample_rate,
+        dtype=int,
     )
     ilines = np.arange(
-        first_iline, first_iline + iline_step * ni - 1e-10, iline_step, dtype=int,
+        first_iline,
+        first_iline + iline_step * ni - 1e-10,
+        iline_step,
+        dtype=int,
     )
     xlines = np.arange(
-        first_xline, first_xline + xline_step * nx - 1e-10, xline_step, dtype=int,
+        first_xline,
+        first_xline + xline_step * nx - 1e-10,
+        xline_step,
+        dtype=int,
     )
 
     if first_offset is None:
diff --git a/segysak/openzgy/__init__.py b/segysak/openzgy/__init__.py
index f02b4e1..5c96561 100644
--- a/segysak/openzgy/__init__.py
+++ b/segysak/openzgy/__init__.py
@@ -1,1 +1,1 @@
-from ._openzgy import zgy_loader, zgy_writer
\ No newline at end of file
+from ._openzgy import zgy_loader, zgy_writer
diff --git a/segysak/segy/__init__.py b/segysak/segy/__init__.py
index 5a0bb68..48178a4 100644
--- a/segysak/segy/__init__.py
+++ b/segysak/segy/__init__.py
@@ -1,4 +1,10 @@
 from ._segy_headers import segy_header_scrape, segy_bin_scrape, segy_header_scan
 from ._segy_text import get_segy_texthead, put_segy_texthead, create_default_texthead
-from ._segy_loader import segy_loader, well_known_byte_locs, segy_converter
+from ._segy_loader import (
+    segy_loader,
+    well_known_byte_locs,
+    segy_converter,
+    segy_freeloader,
+)
 from ._segy_writer import ncdf2segy, segy_writer, output_byte_locs
+from ._segy_freewriter import segy_freewriter
diff --git a/segysak/segy/_segy_freewriter.py b/segysak/segy/_segy_freewriter.py
new file mode 100644
index 0000000..6ca25a2
--- /dev/null
+++ b/segysak/segy/_segy_freewriter.py
@@ -0,0 +1,285 @@
+import importlib
+import xarray as xr
+import segyio
+import numpy as np
+import pandas as pd
+from more_itertools import split_into, split_when
+from itertools import product
+
+try:
+    has_ipywidgets = importlib.util.find_spec("ipywidgets") is not None
+    if has_ipywidgets:
+        from tqdm.autonotebook import tqdm
+    else:
+        from tqdm import tqdm as tqdm
+except ModuleNotFoundError:
+    from tqdm import tqdm as tqdm
+
+from .._accessor import open_seisnc
+from .._keyfield import CoordKeyField
+from ._segy_text import (
+    create_default_texthead,
+    put_segy_texthead,
+    _clean_texthead,
+    trace_header_map_to_text,
+)
+from ._segy_globals import _ISEGY_MEASUREMENT_SYSTEM
+
+TQDM_ARGS = dict(unit=" traces", desc="Writing to SEG-Y")
+
+
+def _chunk_iterator(ds):
+    # get the chunk selection labels
+    if ds.chunks:
+        ranges = {
+            dim: tuple(split_into(ds[dim].values, ds.chunks[dim])) for dim in ds.chunks
+        }
+
+    else:
+        ranges = {dim: [tuple(ds[dim].values)] for dim in ds.dims}
+    lens = {key: len(val) for key, val in ranges.items()}
+    order = tuple(ranges.keys())
+    ranges = tuple(ranges.values())
+    for selection in product(*map(range, tuple(lens.values()))):
+        yield {order[i]: list(ranges[i][n]) for i, n in enumerate(selection)}
+
+
+def _segy_freewriter(
+    seisnc,
+    segyfile,
+    data_array="data",
+    trace_header_map=None,
+    silent=False,
+    use_text=False,
+    dead_trace_key=None,
+    vert_dimension="twt",
+    **dim_kwargs,
+):
+    """"""
+
+    try:
+        coord_scalar = seisnc.coord_scalar
+    except AttributeError:
+        coord_scalar = 0
+        seisnc.attrs["coord_scalar"] = coord_scalar
+    coord_scalar_mult = np.power(abs(coord_scalar), np.sign(coord_scalar) * -1)
+    seisnc.attrs["coord_scalar_mult"] = coord_scalar_mult
+
+    # create empty trace header map if necessary
+    if trace_header_map is None:
+        trace_header_map = dict()
+    else:  # check that values of thm in ds
+        for key in trace_header_map:
+            try:
+                _ = seisnc[key]
+            except KeyError:
+                raise ValueError("keys of trace_header_map must be in seisnc")
+
+    # remove duplicate dim_kwargs from trace_header_map
+    for key in dim_kwargs:
+        trace_header_map[key] = dim_kwargs[key]
+    thm_vars = tuple(trace_header_map.keys())
+    if dead_trace_key:
+        thm_vars += (dead_trace_key,)
+
+    dim_order = tuple(dim_kwargs.keys())
+
+    # transfrom text if requested
+    if use_text:
+        text = {i + 1: line for i, line in enumerate(seisnc.text.split("\n"))}
+        text = _clean_texthead(text, 80)
+    else:
+        text = None
+
+    # create trace headers
+    trace_layout = seisnc[data_array].isel(**{vert_dimension: 0})
+    trace_headers = pd.DataFrame()
+    for var in thm_vars:
+        trace_headers[var] = (
+            seisnc[var]
+            .broadcast_like(trace_layout)
+            .transpose(*dim_order, transpose_coords=True)
+            .values.ravel()
+        )
+
+    # scale coords
+    if CoordKeyField.cdp_x in trace_header_map:
+        trace_headers[CoordKeyField.cdp_x] = (
+            trace_headers[CoordKeyField.cdp_x] * seisnc.coord_scalar_mult
+        )
+    if CoordKeyField.cdp_y in trace_header_map:
+        trace_headers[CoordKeyField.cdp_y] = (
+            trace_headers[CoordKeyField.cdp_y] * seisnc.coord_scalar_mult
+        )
+
+    trace_headers = trace_headers.astype(np.int32)
+    z0 = int(seisnc[vert_dimension].values[0])
+    ntraces = trace_headers.shape[0]
+    ns = seisnc[vert_dimension].values.size
+    msys = _ISEGY_MEASUREMENT_SYSTEM[seisnc.seis.get_measurement_system()]
+    spec = segyio.spec()
+    spec.format = 1
+    spec.iline = 181
+    spec.xline = 185
+    spec.samples = range(ns)
+
+    if dead_trace_key:
+        dead = ~trace_headers[dead_trace_key].astype(bool)
+        spec.tracecount = dead.sum()
+        live_index = trace_headers.where(dead).dropna().index
+        tracen = pd.Series(
+            index=live_index, data=np.arange(live_index.size), name="tracen"
+        )
+    else:
+        spec.tracecount = ntraces
+        tracen = pd.Series(
+            index=trace_headers.index,
+            name="tracen",
+            data=np.arange(trace_headers.index.size),
+        )
+
+    # build trace headers
+    trace_headers = trace_headers.astype(np.int32)
+    trace_headers = trace_headers.join(
+        tracen.astype(np.float32),
+    )
+
+    seisnc["tracen"] = trace_headers.set_index(list(dim_order))["tracen"].to_xarray()
+    trace_headers = trace_headers.dropna(subset=["tracen"])
+    trace_headers["tracen"] = trace_headers["tracen"].astype(np.int32)
+    trace_headers = trace_headers.set_index("tracen")
+    trace_headers = trace_headers[list(trace_header_map.keys())].rename(
+        columns=trace_header_map
+    )
+
+    # header constants
+    for byte, val in zip(
+        [segyio.su.offset, segyio.su.ns, segyio.su.delrt, segyio.su.scalco],
+        [1, ns, z0, int(seisnc.coord_scalar)],
+    ):
+        if byte not in trace_headers.columns:  # don't override user values
+            trace_headers[byte] = val
+
+    # write out to segy
+    with segyio.create(segyfile, spec) as segyf:
+        pbar = tqdm(total=spec.tracecount, disable=silent)  # , **TQDM_ARGS )
+        for chk in _chunk_iterator(seisnc):
+            _ = chk.pop(vert_dimension)
+            # set order
+            data = seisnc.sel(**chk).transpose(*dim_order, vert_dimension)
+            # create missing trace mask and slicing tools
+            mask = ~data.tracen.isnull().values.ravel()
+            trace_locs = data.tracen.values.ravel()[mask].astype(int)
+            trace_locs = tuple(
+                map(
+                    lambda x: slice(x[0], x[-1] + 1),
+                    split_when(trace_locs, lambda x, y: y != x + 1),
+                )
+            )
+            if data[data_array].size > 0:
+                npdata = (
+                    data[data_array].values.reshape(-1, ns)[mask, :].astype(np.float32)
+                )
+                n1 = 0
+                for slc_set in trace_locs:
+                    segyf.header[slc_set] = trace_headers[slc_set].to_dict(
+                        orient="records"
+                    )
+                    n2 = n1 + slc_set.stop - slc_set.start
+                    segyf.trace[slc_set] = npdata[n1:n2]
+                    n1 = n2
+                    pbar.update(slc_set.stop - slc_set.start)
+        pbar.close()
+
+        segyf.bin.update(
+            # tsort=segyio.TraceSortingFormat.INLINE_SORTING,
+            hdt=int(seisnc.sample_rate * 1000),
+            hns=ns,
+            mfeet=msys,
+            jobid=1,
+            lino=1,
+            reno=1,
+            ntrpr=ntraces,
+            nart=ntraces,
+            fold=1,
+        )
+
+    if not text:
+        header_locs_text = trace_header_map_to_text(trace_header_map)
+
+        # create text header
+        overrides = {
+            #     123456789012345678901234567890123456789012345678901234567890123456
+            7: f"Original File: {seisnc.source_file}",
+        }
+
+        l1 = 40 - len(header_locs_text)
+        overrides.update({l1 + i: line for i, line in enumerate(header_locs_text)})
+
+        put_segy_texthead(
+            segyfile,
+            create_default_texthead(overrides),
+        )
+    else:
+        put_segy_texthead(segyfile, text, line_counter=False)
+
+
+def segy_freewriter(
+    seisnc,
+    segyfile,
+    data_array="data",
+    trace_header_map=None,
+    use_text=False,
+    dead_trace_key=None,
+    vert_dimension="twt",
+    chunk_spec=None,
+    silent=False,
+    **dim_kwargs,
+):
+    """Convert siesnc format (NetCDF4) to SEGY.
+
+    Args:
+        seisnc (xarray.Dataset, string): The input SEISNC file either a path or the in memory xarray.Dataset
+        segyfile (string): The output SEGY file
+        data_array (string): The Dataset variable name for the output volume.
+        trace_header_map (dict, optional): Defaults to None. A dictionary of seisnc variables
+            and byte locations. The variable will be written to the trace headers in the
+            assigned byte location. By default CMP=23, cdp_x=181, cdp_y=185, iline=189,
+            xline=193.
+        use_text (book, optional): Use the seisnc text for the EBCIDC output. This text usally comes from
+            the loaded SEG-Y file and may not match the segysak SEG-Y output. Defaults to False and writes
+            the default segysak EBCIDC
+        dead_trace_key (str): The key for the Dataset variable to use for trace output filter.
+        vert_dimension (str): Data dimension to output, defaults to 'twt' or 'depth' whichever
+            is present.
+        chunk_spec (dict, optional): Xarray open_dataset chunking spec.
+        silent (bool, optional): Turn off progress reporting. Defaults to False.
+        dim_kwargs (int): The dimension/byte location pairs to output dimensions to. The number of dim_kwargs should be
+            equal to the number of dimensions on the output data_array. The sort order will be as per the order passed
+            to the function.
+    """
+    kwargs = dict(
+        data_array=data_array,
+        trace_header_map=trace_header_map,
+        silent=silent,
+        use_text=use_text,
+        dead_trace_key=dead_trace_key,
+        vert_dimension=vert_dimension,
+    )
+
+    if isinstance(seisnc, xr.Dataset):
+        _segy_freewriter(
+            seisnc,
+            segyfile,
+            **kwargs,
+            **dim_kwargs,
+        )
+    else:
+        ncfile = seisnc
+        with open_seisnc(ncfile, chunks=chunk_spec) as seisnc:
+            _segy_freewriter(
+                seisnc,
+                segyfile,
+                **kwargs,
+                **dim_kwargs,
+            )
diff --git a/segysak/segy/_segy_globals.py b/segysak/segy/_segy_globals.py
index b9c0020..92d126e 100644
--- a/segysak/segy/_segy_globals.py
+++ b/segysak/segy/_segy_globals.py
@@ -1,4 +1,3 @@
-
 from collections import defaultdict
 
 _SEGY_MEASUREMENT_SYSTEM = defaultdict(lambda: 0)
diff --git a/segysak/segy/_segy_headers.py b/segysak/segy/_segy_headers.py
index 006224f..ecba92e 100644
--- a/segysak/segy/_segy_headers.py
+++ b/segysak/segy/_segy_headers.py
@@ -1,4 +1,5 @@
 import importlib
+import itertools
 import numpy as np
 import pandas as pd
 import segyio
@@ -78,7 +79,14 @@ def segy_bin_scrape(segyfile, **segyio_kwargs):
         return {key: segyf.bin[item] for key, item in bk.items()}
 
 
-def segy_header_scrape(segyfile, partial_scan=None, silent=False, **segyio_kwargs):
+def segy_header_scrape(
+    segyfile,
+    partial_scan=None,
+    silent=False,
+    bytes_filter=None,
+    chunk=100_000,
+    **segyio_kwargs,
+):
     """Scape all data from segy trace headers
 
     Args:
@@ -86,25 +94,52 @@ def segy_header_scrape(segyfile, partial_scan=None, silent=False, **segyio_kwarg
         partial_scan (int): Setting partial scan to a positive int will scan only
             that many traces. Defaults to None.
         silent (bool): Disable progress bar.
+        bytes_filter (list): List of byte locations to load exclusively.
+        chunk (int): Number of traces to read in one go.
 
     Returns:
         pandas.DataFrame: Raw header information in table for scanned traces.
     """
+    assert (chunk > 0) and isinstance(chunk, int)
     header_keys = _active_tracefield_segyio()
-    columns = header_keys.keys()
+    enum_byte_index = {
+        int(byte_loc): i for i, byte_loc in enumerate(header_keys.values())
+    }
+
+    if bytes_filter:
+        for byte_loc in bytes_filter:
+            assert byte_loc in enum_byte_index
+        bytes_filter_index = [enum_byte_index[byte_loc] for byte_loc in bytes_filter]
+        enum_filter = [segyio.TraceField(byte_loc) for byte_loc in bytes_filter]
+    else:
+        bytes_filter_index = list(enum_byte_index.values())
+        enum_filter = [segyio.TraceField(byte_loc) for byte_loc in header_keys.values()]
+
+    columns = [list(header_keys.keys())[i] for i in bytes_filter_index]
     segyio_kwargs["ignore_geometry"] = True
+
     with segyio.open(segyfile, "r", **segyio_kwargs) as segyf:
+        segyf_hgen = segyf.header[:]
         ntraces = segyf.tracecount
         if partial_scan is not None:
-            ntraces = int(partial_scan)
-        slc = slice(0, ntraces, 1)
-        # take headers returned from segyio and create lists for a dataframe
-        hv = map(
-            lambda x: np.array(list(x.values())),
-            tqdm(segyf.header[slc], total=ntraces, disable=silent, **TQDM_ARGS),
-        )
-        head_df = pd.DataFrame(np.vstack(list(hv)), columns=columns)
-        head_df.replace(to_replace=-2147483648, value=np.nan, inplace=True)
+            ntraces = min(ntraces, int(partial_scan))
+
+        head_df = pd.DataFrame(index=pd.Index(range(ntraces)), columns=columns)
+        slc_end = chunk
+        with tqdm(total=ntraces, disable=silent, **TQDM_ARGS) as pbar:
+            while slc_end <= ntraces + chunk - 1:
+                slc = slice(slc_end - chunk, min(slc_end, ntraces), 1)
+                # take headers returned from segyio and create lists for a dataframe
+                head_df.iloc[slc, :] = np.vstack(
+                    [
+                        list(next(segyf_hgen).values())
+                        for _ in range(slc.stop - slc.start)
+                    ]
+                )[:, bytes_filter_index]
+                slc_end += chunk
+                pbar.update(slc.stop - slc.start)
+
+    head_df.replace(to_replace=-2147483648, value=np.nan, inplace=True)
     return head_df
 
 
diff --git a/segysak/segy/_segy_loader.py b/segysak/segy/_segy_loader.py
index 2d570d4..c64f23c 100644
--- a/segysak/segy/_segy_loader.py
+++ b/segysak/segy/_segy_loader.py
@@ -4,6 +4,7 @@
 
 import importlib
 import pathlib
+from numpy.core.fromnumeric import trace
 
 import segyio
 
@@ -38,6 +39,7 @@
     CoordKeyField,
     AttrKeyField,
     VariableKeyField,
+    VerticalKeyDim,
     DimensionKeyField,
 )
 from segysak._seismic_dataset import (
@@ -52,7 +54,12 @@
 
 KNOWN_BYTES = FrozenDict(KNOWN_BYTES)
 
-from ._segy_headers import segy_bin_scrape, segy_header_scrape, what_geometry_am_i
+from ._segy_headers import (
+    segy_bin_scrape,
+    segy_header_scrape,
+    what_geometry_am_i,
+    _active_tracefield_segyio,
+)
 from ._segy_text import get_segy_texthead
 from ._segy_globals import _SEGY_MEASUREMENT_SYSTEM
 
@@ -75,7 +82,7 @@ def _header_to_index_mapping(series):
     as_unique = series.unique()
     sorted_unique = np.sort(as_unique)
     mapper = {val: i for i, val in enumerate(sorted_unique)}
-    return series.replace(mapper)
+    return series.map(mapper.get)
 
 
 def _segy3d_ncdf(
@@ -131,7 +138,9 @@ def _segy3d_ncdf(
         for _, grp in head_df.groupby(contig_dir):
             for trc, val in grp.iterrows():
                 i1, i2 = val[["il_index", "xl_index"]].values.astype(int)
-                seisnc_data[i1, i2, :] = segyf.trace[trc][n0 : ns + 1]
+                seisnc_data[i1, i2, :] = segyf.trace[trc][n0 : ns + 1].astype(
+                    np.float32
+                )
                 percentiles = (
                     np.percentile(segyf.trace[trc][n0 : ns + 1], PERCENTILES)
                     + percentiles
@@ -372,8 +381,8 @@ def _3dsegy_loader(
     xlines = np.sort(xlines)
     iline_index_map = {il: i for i, il in enumerate(ilines)}
     xline_index_map = {xl: i for i, xl in enumerate(xlines)}
-    head_df.loc[:, "il_index"] = head_df[head_loc.iline].replace(iline_index_map)
-    head_df.loc[:, "xl_index"] = head_df[head_loc.xline].replace(xline_index_map)
+    head_df.loc[:, "il_index"] = head_df[head_loc.iline].map(iline_index_map.get)
+    head_df.loc[:, "xl_index"] = head_df[head_loc.xline].map(xline_index_map.get)
 
     # binary header translation
     ns = head_bin["Samples"]
@@ -385,7 +394,7 @@ def _3dsegy_loader(
         offsets = head_df[head_loc.offset].unique()
         offsets = np.sort(offsets)
         offset_index_map = {off: i for i, off in enumerate(offsets)}
-        head_df["off_index"] = head_df[head_loc.offset].replace(offset_index_map)
+        head_df["off_index"] = head_df[head_loc.offset].map(offset_index_map.get)
     else:
         offsets = None
 
@@ -766,12 +775,27 @@ def _loader_converter_header_handling(
     silent=False,
     extra_byte_fields=None,
     head_df=None,
+    optimised_load=False,
     **segyio_kwargs,
 ):
 
+    if optimised_load:
+        scrape_bytes = [
+            val for val in [cdp, iline, xline, cdpx, cdpy, offset] if val is not None
+        ]
+        scrape_bytes += list(extra_byte_fields.values())
+        scrape_bytes += [
+            _active_tracefield_segyio()["DelayRecordingTime"],
+            _active_tracefield_segyio()["SourceGroupScalar"],
+        ]
+    else:
+        scrape_bytes = None
+
     if head_df is None:
         # Start by scraping the headers.
-        head_df = segy_header_scrape(segyfile, silent=silent, **segyio_kwargs)
+        head_df = segy_header_scrape(
+            segyfile, silent=silent, bytes_filter=scrape_bytes, **segyio_kwargs
+        )
 
     head_bin = segy_bin_scrape(segyfile, **segyio_kwargs)
     head_loc = AttrDict(
@@ -885,6 +909,138 @@ def _loader_converter_write_headers(
     return ds
 
 
+def segy_freeloader(
+    segyfile,
+    vert_domain="TWT",
+    data_type="AMP",
+    return_geometry=False,
+    silent=False,
+    extra_byte_fields=None,
+    head_df=None,
+    segyio_kwargs=None,
+    **dim_kwargs,
+):
+    """Freeform loader for SEG-Y data. This loader allows you to load SEG-Y into
+    an xarray.Dataset using an arbitrary number of header locations to create
+    othogonal dimensions. This is an eager loader and will transfer the entire
+    SEG-Y and requested header information to memory.
+
+    From the dimension header locations specified the freeloader will try to
+    create a Dataset where each trace is assigned to a dimension.
+
+    Args:
+        segyfile (string): The SEG-Y file/path.
+        vert_domain (str, optional): One of ('TWT', 'DEPTH'). Defaults to 'TWT'.
+        data_type (str, optional): Defaults to "AMP".
+        return_geometry (bool, optional): If true, just returned the empty
+            dataset based upon the calcuated header geometry. Defaults to False.
+        silent (bool, optional): Turn off progress bars. Defaults to False.
+        extra_byte_fields (dict, optional): Additional header information to
+            load into the Dataset. Defaults to None.
+        head_df (pandas.DataFrame): The DataFrame output from `segy_header_scrape`.
+            This DataFrame can be filtered by the user
+            to load select trace sets. Trace loading is based upon the DataFrame index.
+        segyio_kwargs (dict, optional): Extra keyword arguments for segyio.open
+        **dim_kwargs: Dimension names and byte location pairs.
+    """
+    if segyio_kwargs is None:
+        segyio_kwargs = dict()
+
+    if head_df is None:
+        # Start by scraping the headers.
+        head_df = segy_header_scrape(segyfile, silent=silent, **segyio_kwargs)
+
+    head_bin = segy_bin_scrape(segyfile, **segyio_kwargs)
+
+    # get vertical sample ranges
+    n0 = 0
+    nsamp = head_bin["Samples"]
+    ns0 = head_df.DelayRecordingTime.min()
+
+    # binary header translation
+    nsamp = head_bin["Samples"]
+    sample_rate = head_bin["Interval"] / 1000.0
+    msys = _SEGY_MEASUREMENT_SYSTEM[head_bin["MeasurementSystem"]]
+    vert_samples = np.arange(ns0, ns0 + sample_rate * nsamp, sample_rate, dtype=int)
+
+    # creating dimensions and new dataset
+    dims = dict()
+    dim_index_names = list()
+    dim_fields = list()
+    for dim in dim_kwargs:
+        trace_field = str(segyio.TraceField(dim_kwargs[dim]))
+        if trace_field == "Unknown Enum":
+            raise ValueError(f"{dim}:{dim_kwargs[dim]} was not a valid byte header")
+        dim_fields.append(trace_field)
+        as_unique = head_df[trace_field].unique()
+        dims[dim] = np.sort(as_unique)
+        d_map = {dval: i for i, dval in enumerate(as_unique)}
+        index_name = f"{dim}_index"
+        dim_index_names.append(index_name)
+        head_df.loc[:, index_name] = head_df[trace_field].map(d_map)
+
+    if (
+        head_df[dim_index_names].shape
+        != head_df[dim_index_names].drop_duplicates().shape
+    ):
+        raise ValueError(
+            "The selected dimensions results in multiple traces per "
+            "dimension location, add additional dimensions or use "
+            "trace numbering to load as 2D."
+        )
+
+    builder, domain = _dataset_coordinate_helper(vert_samples, vert_domain, **dims)
+    ds = create_seismic_dataset(**builder)
+
+    # getting attributes
+    text = get_segy_texthead(segyfile, **segyio_kwargs)
+    ds.attrs[AttrKeyField.text] = text
+    ds.attrs[AttrKeyField.source_file] = pathlib.Path(segyfile).name
+    ds.attrs[AttrKeyField.measurement_system] = msys
+    ds.attrs[AttrKeyField.sample_rate] = sample_rate
+
+    # map extra byte fields into ds
+    if extra_byte_fields is not None:
+        to_add = list()
+        for name, byte in extra_byte_fields.items():
+            trace_field = str(segyio.TraceField(byte))
+            if trace_field == "Unknown Enum":
+                raise ValueError(f"{name}:{byte} was not a valid byte header")
+            to_add.append(trace_field)
+        to_add = to_add + dim_fields
+        extras = head_df[to_add].set_index(dim_fields).to_xarray()
+        extras = extras.rename_dims(
+            {b: a for a, b in zip(dim_kwargs, dim_fields) if a != b}
+        )
+        for name, xtr in zip(extra_byte_fields, to_add):
+            ds[name] = extras[xtr]
+
+    if return_geometry:
+        # return geometry -> e.g. don't process segy traces
+        return ds
+
+    segyio_kwargs.update(dict(ignore_geometry=True))
+
+    with segyio.open(segyfile, "r", **segyio_kwargs) as segyf:
+
+        segyf.mmap()
+        shape = [ds.dims[d] for d in dim_kwargs] + [vert_samples.size]
+        volume = np.zeros(shape, dtype=np.float32)
+
+        # this can probably be done as a block - leaving for now just incase sorting becomes an issue
+        indexes = tuple([head_df[idx].values for idx in dim_index_names])
+        volume[indexes] = segyf.trace.raw[:][head_df.index.values]
+
+    percentiles = np.percentile(volume, PERCENTILES)
+    ds[VariableKeyField.data] = (
+        list(dim_kwargs) + [VerticalKeyDim[domain]],
+        volume,
+    )
+    ds.attrs[AttrKeyField.percentiles] = list(percentiles)
+
+    return ds
+
+
 def segy_loader(
     segyfile,
     cdp=None,
@@ -912,7 +1068,7 @@ def segy_loader(
             cdp/iline - CDP or Inline axis
             xline - Xline axis
             twt/depth - The vertical axis
-            d4 - Offset/Angle Axis
+            offset - Offset/Angle Axis
         Coordinates:
             iline - The inline numbering
             xline - The xline numbering
@@ -935,7 +1091,7 @@ def segy_loader(
         segyfile (str): Input segy file path
         cdp (int, optional): The CDP byte location, usually 21.
         iline (int, optional): Inline byte location, usually 189
-        xline (int, optional): Cross-line byte location, usally 193
+        xline (int, optional): Cross-line byte location, usually 193
         cdpx (int, optional): UTMX byte location, usually 181
         cdpy (int, optional): UTMY byte location, usually 185
         offset (int, optional): Offset/angle byte location
@@ -1067,7 +1223,7 @@ def segy_converter(
             cdp/iline - CDP or Inline axis
             xline - Xline axis
             twt/depth - The vertical axis
-            d4 - Offset/Angle Axis
+            offset - Offset/Angle Axis
         Coordinates:
             iline - The inline numbering
             xline - The xline numbering
@@ -1092,7 +1248,7 @@ def segy_converter(
             returned in memory as an xarray.Dataset.
         cdp (int, optional): The CDP byte location, usually 21.
         iline (int, optional): Inline byte location, usually 189
-        xline (int, optional): Cross-line byte location, usally 193
+        xline (int, optional): Cross-line byte location, usually 193
         cdpx (int, optional): UTMX byte location, usually 181
         cdpy (int, optional): UTMY byte location, usually 185
         offset (int, optional): Offset/angle byte location
@@ -1133,9 +1289,10 @@ def segy_converter(
         return_geometry=return_geometry,
         silent=silent,
         extra_byte_fields=extra_byte_fields,
+        optimised_load=True,
         **segyio_kwargs,
     )
-
+    print("header_loaded")
     common_args = (segyfile, head_df, head_bin, head_loc)
 
     common_kwargs = dict(
@@ -1149,6 +1306,8 @@ def segy_converter(
 
     # 3d data needs iline and xline
     if iline is not None and xline is not None:
+
+        print("is_3d")
         ds = _3dsegy_loader(
             *common_args,
             **common_kwargs,
diff --git a/segysak/segy/_segy_text.py b/segysak/segy/_segy_text.py
index a29384c..3eb64a4 100644
--- a/segysak/segy/_segy_text.py
+++ b/segysak/segy/_segy_text.py
@@ -67,9 +67,7 @@ def get_segy_texthead(segyfile, ext_headers=False, no_richstr=False, **segyio_kw
         return _upgrade_txt_richstr(text)
 
 
-def put_segy_texthead(
-    segyfile, ebcidc, ext_headers=False, line_counter=True, **segyio_kwargs
-):
+def put_segy_texthead(segyfile, ebcidc, line_counter=True, **segyio_kwargs):
 
     header = ""
     if isinstance(ebcidc, dict):
@@ -106,7 +104,7 @@ def put_segy_texthead(
             warn("String EBCIDC is too large - truncating", UserWarning)
         header = bytes(ebcidc[:3200], "utf8")
     else:
-        raise ValueError("Unknown ebcidc type")
+        raise ValueError("Unknown ebcidc type expect bytes, str or dict")
 
     segyio_kwargs["ignore_geometry"] = True
     with segyio.open(segyfile, "r+", **segyio_kwargs) as segyf:
@@ -182,12 +180,34 @@ def create_default_texthead(override=None):
         4: f"DATA FORMAT: SEG-Y;  DATE: {today} {time}",
         5: "DATA DESCRIPTION: SEG-Y format data output from segysak using segyio",
         6: "",
-        35: "*** BYTE LOCATION OF KEY HEADERS ***",
-        36: "CMP UTM-X 181-184, ALL COORDS X100, CMP UTM-Y 185-188",
-        37: "INLINE 189-193, XLINE 194-198, ",
         40: "END TEXTUAL HEADER",
     }
     if override is not None:
         for key, line in override.items():
             text_dict[key] = line
     return _clean_texthead(text_dict)
+
+
+def trace_header_map_to_text(trace_header_map):
+    """Convert a trace header map to lines of text to say where header info was
+    put in a file.
+
+    Args:
+        trace_header_map (dict): Header byte locations.
+
+    Returns:
+        list: List of lines to output.
+    """
+    header_locs_text = [f"{key}: {value}," for key, value in trace_header_map.items()]
+
+    lines = ["*** BYTE LOCATION OF KEY HEADERS ***"]
+    line = ""
+    for byte_loc in header_locs_text:
+        if len(line + byte_loc) < 79:
+            line += byte_loc
+            line += " "
+        else:
+            lines.append(line)
+            line = byte_loc
+    lines.append(line)
+    return lines
diff --git a/segysak/segy/_segy_writer.py b/segysak/segy/_segy_writer.py
index 6c96551..80ced66 100644
--- a/segysak/segy/_segy_writer.py
+++ b/segysak/segy/_segy_writer.py
@@ -261,7 +261,7 @@ def _ncdf2segy_3d(
     silent=False,
     **thm_args,
 ):
-    """Convert etlpy siesnc format (NetCDF4) to SEGY.
+    """Convert siesnc format (NetCDF4) to SEGY.
 
     Args:
         ds (xarray.Dataset): The input SEISNC dataset
@@ -298,7 +298,7 @@ def _ncdf2segy_3d(
     # the file, i.e. its inline numbers, crossline numbers, etc. You can also add
     # more structural information, but offsets etc. have sensible defautls. This is
     # the absolute minimal specification for a N-by-M volume
-    spec.sorting = 1
+    spec.sorting = 2
     spec.format = 1
     spec.iline = iline
     spec.xline = xline
@@ -332,19 +332,6 @@ def _ncdf2segy_3d(
                 pbar.update(nj)
         pbar.close()
 
-        segyf.bin.update(
-            tsort=segyio.TraceSortingFormat.INLINE_SORTING,
-            hdt=int(ds.sample_rate * 1000),
-            hns=nk,
-            mfeet=msys,
-            jobid=1,
-            lino=1,
-            reno=1,
-            ntrpr=ni * nj,
-            nart=ni * nj,
-            fold=1,
-        )
-
     if not text:
         # create text header
         overrides = {
@@ -378,7 +365,7 @@ def _ncdf2segy_3dgath(
     silent=False,
     **thm_args,
 ):
-    """Convert etlpy siesnc format (NetCDF4) to SEGY.
+    """Convert siesnc format (NetCDF4) to SEGY.
 
     Args:
         ds (xarray.Dataset): The input SEISNC dataset
@@ -418,7 +405,7 @@ def _ncdf2segy_3dgath(
     # the file, i.e. its inline numbers, crossline numbers, etc. You can also add
     # more structural information, but offsets etc. have sensible defautls. This is
     # the absolute minimal specification for a N-by-M volume
-    spec.format = 1
+    spec.format = 2
     spec.iline = iline
     spec.xline = xline
     spec.samples = ds[dimension].astype(int).values
@@ -495,7 +482,7 @@ def _ncdf2segy_2d(
     silent=False,
     **thm_args,
 ):
-    """Convert etlpy siesnc format (NetCDF4) to SEGY.
+    """Convert siesnc format (NetCDF4) to SEGY.
 
     Args:
         ds (xarray.Dataset): The input SEISNC dataset
@@ -595,7 +582,7 @@ def _ncdf2segy_2d_gath(
     silent=False,
     **thm_args,
 ):
-    """Convert etlpy siesnc format (NetCDF4) to SEGY.
+    """Convert siesnc format (NetCDF4) to SEGY.
 
     Args:
         ds (xarray.Dataset): The input SEISNC dataset
@@ -724,6 +711,7 @@ def _segy_writer_input_handler(
         _ncdf2segy_2d(ds, segyfile, **common_kwargs, **thm)
     elif ds.seis.is_2dgath():
         thm = output_byte_locs("standard_2d_gath")
+        thm.update(trace_header_map)
         _ncdf2segy_2d_gath(ds, segyfile, **common_kwargs, **thm)
     elif ds.seis.is_3d():
         thm = output_byte_locs("standard_3d")
diff --git a/segysak/seisnc.py b/segysak/seisnc.py
index 80d5976..30be9e0 100644
--- a/segysak/seisnc.py
+++ b/segysak/seisnc.py
@@ -115,9 +115,30 @@ def set_seisnc_dims(
 
     # update seisnc on disk
     with netCDF4.Dataset(str(ncfile), mode="a") as seisnc:
-        seisnc.createVariable("CDP_X", "f4", ("i", "j",))
-        seisnc.createVariable("CDP_Y", "f4", ("i", "j",))
-        seisnc.createVariable("CDP_TRACE", "i4", ("i", "j",))
+        seisnc.createVariable(
+            "CDP_X",
+            "f4",
+            (
+                "i",
+                "j",
+            ),
+        )
+        seisnc.createVariable(
+            "CDP_Y",
+            "f4",
+            (
+                "i",
+                "j",
+            ),
+        )
+        seisnc.createVariable(
+            "CDP_TRACE",
+            "i4",
+            (
+                "i",
+                "j",
+            ),
+        )
         seisnc.createVariable("data", "f4", ("i", "j", "k"))
 
         seisnc.createVariable("vert", "f4", ("k"))
diff --git a/segysak/tools.py b/segysak/tools.py
index b209501..1124fd8 100644
--- a/segysak/tools.py
+++ b/segysak/tools.py
@@ -173,4 +173,3 @@ def halfsample(arr):
 def plane(xy, a, b, c):
     x, y = xy
     return a * x + b * y + c
-
diff --git a/setup.py b/setup.py
index 45c3feb..1aaeb06 100644
--- a/setup.py
+++ b/setup.py
@@ -22,6 +22,7 @@
     "h5netcdf",
     "setuptools_scm",
     "attrdict",
+    "more-itertools",
 ]
 
 notebook_deps = [
@@ -33,7 +34,14 @@
     "pyvista",
 ]
 
-testing_deps = ["flake8", "pytest", "hypothesis", "affine", "pytest_cases"]
+testing_deps = [
+    "flake8",
+    "pytest",
+    "hypothesis",
+    "affine",
+    "pytest_cases",
+    "subsurface",
+]
 
 extras_require = {
     "notebook": notebook_deps,
@@ -49,6 +57,7 @@
         "sphinx-copybutton",
         "rtds-action",
         "jupytext",
+        "pyvista",
     ]
     + notebook_deps,
 }


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index c74d072..5bff47c 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -4,6 +4,7 @@
 from .fixtures_segy import *
 from .fixtures_zgy import *
 from .fixtures_seisnc import *
+from .fixtures_data import *
 
 collect_ignore = ["test_seismic_zgy.py"]
 
diff --git a/tests/fixtures_data.py b/tests/fixtures_data.py
new file mode 100644
index 0000000..296761b
--- /dev/null
+++ b/tests/fixtures_data.py
@@ -0,0 +1,65 @@
+import pathlib
+import pytest
+from segysak.segy import segy_loader
+
+TESTS_PATH = pathlib.Path(__file__).parent.absolute()
+TEST_DATA_SEGYIO = TESTS_PATH / "test-data-segyio"
+TEST_DATA_SEGYSAK = TESTS_PATH / "test-data-segysak"
+TEST_DATA_VOLVE = TESTS_PATH / "test-data-volve"
+
+
+@pytest.fixture(scope="session")
+def f3_dataset():
+    ds = segy_loader(
+        TEST_DATA_SEGYIO / "f3.sgy",
+        **{"iline": 189, "xline": 193, "cdpx": 181, "cdpy": 185},
+        silent=True
+    )
+
+    ds["extra"] = ds.cdp_x * 0.0
+    return ds
+
+
+@pytest.fixture(scope="session")
+def f3_horizon_exact(f3_dataset):
+    df = f3_dataset.drop_vars("data").to_dataframe().reset_index()
+    df["horizon"] = 50.0
+    return df
+
+
+@pytest.fixture(scope="session")
+def f3_horizon_shift(f3_dataset):
+    df = f3_dataset.drop_vars("data").to_dataframe().reset_index()
+    df["horizon"] = 50.0
+    df["cdp_x"] = df["cdp_x"] + 3.5
+    df["cdp_y"] = df["cdp_y"] + 3.5
+    return df
+
+
+@pytest.fixture(scope="session")
+def geometry_dataset():
+    ds = segy_loader(TEST_DATA_SEGYSAK / "test-geometry.sgy", silent=True)
+    return ds
+
+
+@pytest.fixture(scope="session")
+def f3_withdead_dataset():
+    ds = segy_loader(
+        TEST_DATA_SEGYSAK / "f3-withdead.sgy",
+        **{"iline": 189, "xline": 193, "cdpx": 181, "cdpy": 185},
+        silent=True
+    )
+    return ds
+
+
+@pytest.fixture(scope="session")
+def volve_2d_dataset():
+    ds = segy_loader(
+        TEST_DATA_VOLVE / "volve10r12-full-z-il10117sub.sgy",
+        silent=True,
+        cdp=21,
+        cdpx=73,
+        cdpy=77,
+        vert_domain="DEPTH",
+    )
+    return ds
diff --git a/tests/fixtures_segy.py b/tests/fixtures_segy.py
index aa37a33..8fc09fd 100644
--- a/tests/fixtures_segy.py
+++ b/tests/fixtures_segy.py
@@ -15,8 +15,9 @@
 TEST_SEGY_REG = "test_reg.segy"
 TEST_SEGY_SKEW = "test_skew.segy"
 
-TEST_DATA_SEGYIO = pathlib.Path(__file__).parent.absolute() / "test-data-segyio"
-TEST_DATA_SEGYIO = pathlib.Path(TEST_DATA_SEGYIO)
+TEST_DATA_PATH = pathlib.Path(__file__).parent.absolute()
+TEST_DATA_SEGYIO = TEST_DATA_PATH / "test-data-segyio"
+TEST_DATA_SEGYSAK = TEST_DATA_PATH / "test-data-segysak"
 
 TEST_FILES_SEGYIO_3D = [
     ("acute-small.sgy", {"iline": 189, "xline": 193}),
@@ -187,3 +188,10 @@ def segyio3dps_test_files(request):
 def segyio_nohead_test_files(request):
     file, segyio_kwargs = request.param
     return (file, segyio_kwargs)
+
+
+@pytest.fixture(
+    scope="session",
+)
+def segy_with_nonascii_text():
+    return TEST_DATA_SEGYSAK / "f3-ebcidctest.sgy"
diff --git a/tests/test-data-segysak/f3-ebcidctest.sgy b/tests/test-data-segysak/f3-ebcidctest.sgy
new file mode 100644
index 0000000..c4284d2
Binary files /dev/null and b/tests/test-data-segysak/f3-ebcidctest.sgy differ
diff --git a/tests/test-data-segysak/f3-withdead.sgy b/tests/test-data-segysak/f3-withdead.sgy
new file mode 100644
index 0000000..6f192c0
Binary files /dev/null and b/tests/test-data-segysak/f3-withdead.sgy differ
diff --git a/tests/test-data-segysak/test-geometry.sgy b/tests/test-data-segysak/test-geometry.sgy
new file mode 100644
index 0000000..7d82935
Binary files /dev/null and b/tests/test-data-segysak/test-geometry.sgy differ
diff --git a/tests/test-data-volve/volve10r12-full-z-il10117sub.sgy b/tests/test-data-volve/volve10r12-full-z-il10117sub.sgy
new file mode 100644
index 0000000..8b0b7da
Binary files /dev/null and b/tests/test-data-volve/volve10r12-full-z-il10117sub.sgy differ
diff --git a/tests/test_accessor.py b/tests/test_accessor.py
index 2a1e194..696ba42 100644
--- a/tests/test_accessor.py
+++ b/tests/test_accessor.py
@@ -1,5 +1,9 @@
 import pytest
+
+import sys
+
 import numpy as np
+from scipy import interpolate
 import xarray as xr
 
 from hypothesis import given, assume, settings
@@ -8,11 +12,168 @@
 
 from affine import Affine
 
+if sys.version_info >= (3, 8):
+    from subsurface import StructuredData
+
 from segysak import create3d_dataset
 from segysak._keyfield import DimensionKeyField
 
+from segysak._accessor import coordinate_df
+
+
+class TestSeisIO:
+    def test_to_subsurface_3d(self, empty3d):
+        if sys.version_info >= (3, 8):
+            assert isinstance(empty3d.seisio.to_subsurface(), StructuredData)
+        else:
+            print("Python <=3.7 not support by subsurface")
+            assert True
+
+    def test_to_subsurface_3dgath(self, empty3d_gath):
+        if sys.version_info >= (3, 8):
+            with pytest.raises(NotImplementedError):
+                empty3d_gath.seisio.to_subsurface()
+        else:
+            print("Python <=3.7 not support by subsurface")
+            assert True
+
+    def test_to_subsurface_2dgath(self, empty2d_gath):
+        if sys.version_info >= (3, 8):
+            with pytest.raises(NotImplementedError):
+                empty2d_gath.seisio.to_subsurface()
+        else:
+            print("Python <=3.7 not support by subsurface")
+            assert True
+
+
+@pytest.mark.parametrize(
+    "coord,extras,linear_fillna,expected_shape",
+    [
+        (True, None, True, (414, 4)),
+        (False, None, True, (414, 4)),
+        (False, None, False, (414, 4)),
+        (
+            True,
+            [
+                "extra",
+            ],
+            True,
+            (414, 5),
+        ),
+    ],
+)
+def test_coordinate_df(f3_dataset, coord, extras, linear_fillna, expected_shape):
+    print(f3_dataset)
+    df = coordinate_df(
+        f3_dataset, coord=coord, extras=extras, linear_fillna=linear_fillna
+    )
+    assert df.shape == expected_shape
+
+
+def test_humanbytes(f3_dataset):
+    assert isinstance(f3_dataset.seis.humanbytes, str)
+
+
+def test_get_measurement_system(f3_dataset):
+    assert f3_dataset.seis.get_measurement_system() == "m"
+    local = f3_dataset.copy(deep=True)
+    local.attrs = {}
+    assert local.seis.get_measurement_system() is None
+
+
+def test_surface_from_points_nparray(f3_dataset, f3_horizon_exact):
+    interpolated = f3_dataset.seis.surface_from_points(
+        f3_horizon_exact[["cdp_x", "cdp_y"]].values,
+        f3_horizon_exact["horizon"].values,
+    )
+    assert interpolated.data.mean() == 50.0
+
+
+def test_surface_from_points_df(f3_dataset, f3_horizon_shift):
+    interpolated = f3_dataset.seis.surface_from_points(
+        f3_horizon_shift, "horizon", right=("cdp_x", "cdp_y")
+    )
+    assert interpolated.horizon.mean() == 50.0
+
+
+def test_subsample_dims(f3_dataset):
+    ss = f3_dataset.seis.subsample_dims(twt=2, xline=1)
+    assert ss["twt"].size == (f3_dataset.twt.values.size * 2 ** 2 - 3)
+    assert ss["xline"].size == (f3_dataset.xline.values.size * 2 - 1)
+
+
+def test_fill_cdpna(f3_dataset):
+    ds = f3_dataset.copy(deep=True)
+    ds.cdp_x[:, 5:10] = np.nan
+    ds.cdp_y[:, 5:10] = np.nan
+    ds.seis.fill_cdpna()
+    assert np.allclose(f3_dataset.cdp_x.values, ds.cdp_x.values, atol=0.01)
+    assert np.allclose(f3_dataset.cdp_y.values, ds.cdp_y.values, atol=0.01)
+
+
+def test_get_affine_transform(geometry_dataset):
+    at = geometry_dataset.seis.get_affine_transform()
+    df = geometry_dataset.drop_vars("data").to_dataframe().reset_index()
+    assert np.allclose(
+        at.transform(df[["iline", "xline"]].values),
+        df[["cdp_x", "cdp_y"]].values,
+        atol=1,
+    )
+
+
+def test_calc_corner_points(f3_dataset):
+    f3_dataset.seis.calc_corner_points()
+    print(f3_dataset.attrs["corner_points"])
+    print(f3_dataset.attrs["corner_points_xy"])
+
+    assert f3_dataset.attrs["corner_points"] == (
+        (111, 875),
+        (111, 892),
+        (133, 892),
+        (133, 875),
+    )
+    assert np.allclose(
+        np.array(f3_dataset.attrs["corner_points_xy"]),
+        np.array(
+            (
+                (620197.2, 6074233.0),
+                (620622.1, 6074245.0),
+                (620606.7, 6074794.5),
+                (620181.94, 6074782.5),
+            )
+        ),
+    )
+
+
+def test_calc_corner_points2d(volve_2d_dataset):
+    volve_2d_dataset.seis.calc_corner_points()
+    assert volve_2d_dataset.attrs["corner_points"] == (1, 202)
+    assert np.allclose(
+        np.array(volve_2d_dataset.attrs["corner_points_xy"]),
+        np.array(
+            (
+                (436482.16, 6477774.5),
+                (436482.16, 6477774.5),
+                (434044.3, 6477777.0),
+                (434044.3, 6477777.0),
+            )
+        ),
+    )
+
+
+def test_get_dead_trace_map(f3_withdead_dataset):
+    dead = f3_withdead_dataset.seis.get_dead_trace_map(zeros_as_nan=True)
+    assert isinstance(dead, xr.DataArray)
+    assert dead.sum().values == 175
+
+    dead = f3_withdead_dataset.seis.get_dead_trace_map(
+        scan=[30, 31, 32], zeros_as_nan=True
+    )
+    assert isinstance(dead, xr.DataArray)
+    assert dead.sum().values == 175
+
 
-def TestIsFunctions():
+class TestIsFunctions:
     def test_is_2d(self, empty2d):
         assert empty2d.seis.is_2d()
         assert not empty2d.seis.is_3d()
diff --git a/tests/test_core.py b/tests/test_core.py
index f05c717..490c765 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -45,4 +45,3 @@ def test_repr(self):
     def test_fromkeys(self):
         fdict = FrozenDict().fromkeys(["a"], [1])
         assert isinstance(fdict, FrozenDict)
-
diff --git a/tests/test_richstring.py b/tests/test_richstring.py
new file mode 100644
index 0000000..3a14104
--- /dev/null
+++ b/tests/test_richstring.py
@@ -0,0 +1,11 @@
+import pytest
+
+from segysak._richstr import _upgrade_txt_richstr
+
+
+def test_richstring(volve_2d_dataset):
+    richstr = _upgrade_txt_richstr(volve_2d_dataset.text)
+    html = richstr._repr_html_()
+    print(richstr)
+    print(html)
+    assert True
diff --git a/tests/test_segy_segy_freewriter.py b/tests/test_segy_segy_freewriter.py
new file mode 100644
index 0000000..8e53504
--- /dev/null
+++ b/tests/test_segy_segy_freewriter.py
@@ -0,0 +1,24 @@
+import pytest
+
+from segysak.segy._segy_freewriter import _chunk_iterator, segy_freewriter
+
+
+def test_chunk_iterator(f3_dataset):
+    ds = f3_dataset.chunk({"iline": 2, "xline": 2})
+    chunks = [chk for chk in _chunk_iterator(ds)]
+    assert len(chunks) == 12 * 9
+
+    chunks = [chk for chk in _chunk_iterator(f3_dataset)]
+    assert len(chunks) == 1
+
+
+def test_segy_freewriter_ds(f3_dataset, temp_dir):
+    sgyf = temp_dir / "f3_test.sgy"
+    segy_freewriter(
+        f3_dataset,
+        str(sgyf),
+        trace_header_map=dict(cdp_x=181, cdp_y=185),
+        iline=189,
+        xline=193,
+    )
+    assert sgyf.exists()
diff --git a/tests/test_seismic_dataset.py b/tests/test_seismic_dataset.py
index 300e9ad..822a95a 100644
--- a/tests/test_seismic_dataset.py
+++ b/tests/test_seismic_dataset.py
@@ -110,6 +110,24 @@ def test_create_2D_seismic_dataset_with_multiple_dimensions(self, d):
         )
         assert len(dataset.dims) == d + 2
 
+    def test_mutally_not_allowed_arguments(self):
+        with pytest.raises(ValueError):
+            ds = create_seismic_dataset(cdp=100, iline=100, xline=100)
+
+        with pytest.raises(ValueError):
+            ds = create_seismic_dataset(cdp=100, iline=100)
+
+        with pytest.raises(ValueError):
+            ds = create_seismic_dataset(cdp=100, xline=100)
+
+    def test_mutually_required_arguments(self):
+
+        with pytest.raises(ValueError):
+            ds = create_seismic_dataset(iline=100)
+
+        with pytest.raises(ValueError):
+            ds = create_seismic_dataset(xline=100)
+
 
 class TestCreate2DDataset:
     """
@@ -128,6 +146,18 @@ def test_create_2D_dataset_custom_cdp(self, t, f, s):
         dataset = create2d_dataset(dims=(t, 100), first_cdp=f, cdp_step=s)
         assert dataset.cdp.data.max() == f + s * t - s
 
+    @given(integers(1, 10000), integers(0, 100), integers(1, 100), integers(0, 50))
+    def test_create_2D_dataset_wfirstoffset(self, s, f, r, o):
+        dataset = create2d_dataset(
+            dims=(100, s, 5),
+            first_cdp=f,
+            cdp_step=s,
+            sample_rate=r,
+            first_offset=o,
+            offset_step=10,
+        )
+        assert dataset.offset.data.max() == 4 * 10 + o
+
 
 class TestCreate3DDataset:
     """
diff --git a/tests/test_seismic_segy.py b/tests/test_seismic_segy.py
index f2adb0b..535e0cf 100644
--- a/tests/test_seismic_segy.py
+++ b/tests/test_seismic_segy.py
@@ -63,13 +63,6 @@ def cropping_limits():
     return crop, cropz
 
 
-def test_get_segy_texthead(temp_dir, temp_segy):
-    ebcidc = get_segy_texthead(temp_segy)
-    assert isinstance(ebcidc, str)
-    assert len(ebcidc) < 3200
-    assert ebcidc[:3] == "C01"
-
-
 def test_well_known_byte_locs_fail():
     with pytest.raises(ValueError):
         well_known_byte_locs("gibberish")
@@ -96,6 +89,9 @@ def test_segy_header_scan(temp_dir, temp_segy):
     assert scanned == TEST_SEGY_SIZE
     # assert isinstance(head, dict) now a dataframe
 
+    with pytest.raises(ValueError):
+        header = segy_header_scan("", max_traces_scan=23.9, silent=True)
+
 
 def test_segy_header_scan_all(temp_dir, temp_segy):
     _ = segy_header_scan(temp_segy, max_traces_scan=0, silent=True)
diff --git a/tests/test_seismic_segy_text.py b/tests/test_seismic_segy_text.py
new file mode 100644
index 0000000..932da4c
--- /dev/null
+++ b/tests/test_seismic_segy_text.py
@@ -0,0 +1,119 @@
+from numpy import isin
+import pytest
+import shutil
+
+from segysak.segy._segy_text import (
+    put_segy_texthead,
+    get_segy_texthead,
+    _isascii,
+    _clean_texthead,
+    create_default_texthead,
+)
+
+
+@pytest.fixture
+def fake_header():
+    fake_header = ""
+    for i in range(40):
+        fake_header += f"C{i:02d} A SEGY HEADER" + " " * 63
+    return fake_header
+
+
+@pytest.fixture
+def default_header():
+    head_dict = create_default_texthead(
+        {20: "Hello World!", 21: "A Very Long Line Repeated" * 10}
+    )
+    return head_dict
+
+
+@pytest.fixture
+def ascii_bytes_header(fake_header):
+    return bytes(fake_header, "ascii")
+
+
+@pytest.fixture
+def segy_ascii(temp_dir, ascii_bytes_header, segy_with_nonascii_text):
+    name = temp_dir / "ascii_segy.sgy"
+    shutil.copy(segy_with_nonascii_text, name)
+    with open(name, mode="wb") as segyf:
+        segyf.seek(0, 0)
+        segyf.write(ascii_bytes_header)
+
+    return name
+
+
+def test_isascii(
+    ascii_bytes_header,
+):
+    assert _isascii(ascii_bytes_header)
+
+
+def test_get_segy_texthead(temp_dir, temp_segy):
+    ebcidc = get_segy_texthead(temp_segy)
+    ebcidc = get_segy_texthead(temp_segy, no_richstr=True)
+
+    assert isinstance(ebcidc, str)
+    assert len(ebcidc) < 3200
+    assert ebcidc[:3] == "C01"
+
+
+def test_get_segy_texthead_odd(temp_dir, segy_with_nonascii_text):
+    ebcidc = get_segy_texthead(segy_with_nonascii_text)
+    print(ebcidc)
+    assert isinstance(ebcidc, str)
+
+
+def test_get_segy_texthead_ascii(segy_ascii):
+    ebcidc = get_segy_texthead(segy_ascii)
+    print(ebcidc)
+    assert isinstance(ebcidc, str)
+
+
+# TODO: NEED AN EXAMPLE WITH MULTIPLE EXTENDED EBCIDC HEADERS
+
+
+def test_clean_texthead(default_header):
+    cleaned = _clean_texthead(default_header)
+    for key, val in cleaned.items():
+        assert isinstance(key, int)
+        assert len(val) <= 80
+
+
+def test_put_segy_texthead(temp_segy, default_header):
+    put_segy_texthead(temp_segy, default_header)
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "Hello World" in ebcidc
+    put_segy_texthead(temp_segy, default_header, line_counter=False)
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "Hello World" in ebcidc
+
+    default_header.update({34.5: "Bad Key", 30: "Long Line" * 20})
+    put_segy_texthead(temp_segy, default_header)
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "Hello World" in ebcidc
+
+
+def test_put_segy_texthead_bytes(temp_segy, ascii_bytes_header):
+    put_segy_texthead(temp_segy, ascii_bytes_header)
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "A SEGY HEADER" in ebcidc
+
+    put_segy_texthead(temp_segy, ascii_bytes_header * 10)
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "A SEGY HEADER" in ebcidc
+
+
+def test_put_segy_texthead_str(temp_segy, ascii_bytes_header):
+    put_segy_texthead(temp_segy, ascii_bytes_header.decode("ascii"))
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "A SEGY HEADER" in ebcidc
+
+    put_segy_texthead(temp_segy, ascii_bytes_header.decode("ascii") * 2)
+    ebcidc = get_segy_texthead(temp_segy)
+    assert "A SEGY HEADER" in ebcidc
+
+
+def test_put_segy_texthead_unknown(temp_segy):
+    with pytest.raises(ValueError):
+        _ = put_segy_texthead(temp_segy, None)
diff --git a/tests/test_tools.py b/tests/test_tools.py
index 0950a23..6724bd0 100644
--- a/tests/test_tools.py
+++ b/tests/test_tools.py
@@ -14,6 +14,21 @@ def test_tools_check_crop():
     assert [5, 10, 25, 28] == tools.check_crop(IN_CROP, IN_LIMITS)
 
 
+def test_tools_check_crop_errors():
+    with pytest.raises(ValueError):
+        tools.check_crop([100, 0, 0, 100], [0, 100, 0, 100])
+    with pytest.raises(ValueError):
+        tools.check_crop([0, 100, 0, 100], [100, 0, 0, 100])
+    with pytest.raises(ValueError):
+        tools.check_crop(
+            [0, 100, 100, 0],
+            [0, 100, 0, 100],
+        )
+
+    with pytest.raises(ValueError):
+        tools.check_crop([10, 100, 10, 100], [0, 9, 0, 9])
+
+
 @given(floats(-10_000, 10_000), floats(-10_000, 10_000), integers(2, 1000))
 def test_halfsample(a, b, n):
     test_array = np.linspace(a, b, n)
@@ -21,6 +36,23 @@ def test_halfsample(a, b, n):
     assert np.allclose(tools.halfsample(test_array), truth_array)
 
 
+def test_check_zcrop():
+    assert [0, 1000] == tools.check_zcrop([0, 1000], [0, 2000])
+
+
+def test_fix_bad_chars():
+    text = "abcdefg�"
+    assert tools.fix_bad_chars(text) == "abcdefg#"
+
+
+def test_get_datetime():
+    assert isinstance(tools._get_datetime(), tuple)
+
+
+def test_get_userid():
+    assert isinstance(tools._get_userid(), str)
+
+
 if __name__ == "__main__":
     c = tools.check_crop(IN_CROP, IN_LIMITS)
     print(c)