Repository: fatiando/verde (python)

Issue:
Function to read in a grid in Surfer format
Fatiando has the [`load_surfer`](http://www.fatiando.org/dev/api/fatiando.datasets.load_surfer.html#fatiando.datasets.load_surfer) function for reading in a grid in [Golden Software's Surfer](http://www.goldensoftware.com/products/surfer) ascii format. We should port that function here but make it return a xarray.Dataset or DataArray instead.

Pull request:
Surfer reader feature
ported the surfer reader from fatiando and changed the output to a xarray Data Array. Ran tests and formatted locally. 

I ported the surfer reader feature from fatiando and changed the output from a dict to an xarray data array. Will you look through it to make sure it's what you had in mind for this feature. Do you want me to add examples to the docstring as well? Right now it just has the docstring from fatiando modified as appropriate.

Fixes #38

**Reminders**

- [x] Run `make format` and `make check` to make sure the code follows the style guide.
- [x] Add tests for new features or tests that would have caught the bug that you're fixing.
- [x] Add new public functions/methods/classes to `doc/api/index.rst` and `verde/__init__.py`.
- [x] Write detailed docstrings for all functions/classes/methods.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.


Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index b2a99d2..7ef2a85 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -73,6 +73,15 @@ Utilities
     grid_to_table
     median_distance
 
+Input/Output
+------------
+
+.. autosummary::
+   :toctree: generated/
+
+    load_surfer
+
+
 .. automodule:: verde.datasets
 
 .. currentmodule:: verde
diff --git a/verde/__init__.py b/verde/__init__.py
index 2ff1b89..6092c57 100644
--- a/verde/__init__.py
+++ b/verde/__init__.py
@@ -14,6 +14,7 @@ from .coordinates import (
 )
 from .mask import distance_mask
 from .utils import variance_to_weights, maxabs, grid_to_table
+from .io import load_surfer
 from .distances import median_distance
 from .blockreduce import BlockReduce, BlockMean
 from .scipygridder import ScipyGridder
diff --git a/verde/io.py b/verde/io.py
new file mode 100644
index 0000000..aae7553
--- /dev/null
+++ b/verde/io.py
@@ -0,0 +1,125 @@
+"""
+Functions for input and output of grids in less common formats.
+"""
+import numpy as np
+import xarray as xr
+
+
+def load_surfer(fname, dtype="float64"):
+    """
+    Read data from a Surfer ASCII grid file.
+
+    This function reads a Surfer grid file, masks any NaNs in the data,
+    and outputs a :class:`xarray.DataArray` that contains easting and northing
+    coordinates, data values, and associated file metadata.
+
+    Surfer is a contouring, griding and surface mapping software
+    from GoldenSoftware. The names and logos for Surfer and Golden
+    Software are registered trademarks of Golden Software, Inc.
+
+    http://www.goldensoftware.com/products/surfer
+
+    Parameters
+    ----------
+    fname : str or file-like object
+        Name or path of the Surfer grid file or an open file (or file-like) object.
+    dtype : str or numpy.dtype
+        The type of variable used for the data. Default is ``float64``. Use
+        ``float32`` if the data are large and precision is not an issue.
+
+    Returns
+    ----------
+    data : :class:`xarray.DataArray`
+        A 2D grid with the data.
+
+    """
+    # Surfer ASCII grid structure
+    # DSAA            Surfer ASCII GRD ID
+    # nCols nRows     number of columns and rows
+    # xMin xMax       X min max
+    # yMin yMax       Y min max
+    # zMin zMax       Z min max
+    # z11 z21 z31 ... List of Z values
+
+    # Only open a file if given a path instead of a file-like object
+    ispath = not hasattr(fname, "readline")
+    if ispath:
+        input_file = open(fname, "r")
+    else:
+        input_file = fname
+    try:
+        grid_id, shape, region, data_range = _read_surfer_header(input_file)
+        field = np.loadtxt(input_file, dtype=dtype)
+        nans = field >= 1.70141e38
+        if np.any(nans):
+            field = np.ma.masked_where(nans, field)
+        _check_surfer_integrity(field, shape, data_range)
+        attrs = {"gridID": grid_id}
+        if ispath:
+            attrs["file"] = fname
+        dims = ("northing", "easting")
+        coords = {
+            "northing": np.linspace(*region[2:], shape[0]),
+            "easting": np.linspace(*region[:2], shape[1]),
+        }
+        data = xr.DataArray(field, coords=coords, dims=dims, attrs=attrs)
+    finally:
+        if ispath:
+            input_file.close()
+    return data
+
+
+def _read_surfer_header(input_file):
+    """
+    Parse the header record of the grid file.
+
+    The header contains information on the grid shape, region, and the minimum and
+    maximum data values.
+
+    Parameters
+    ----------
+    input_file : file-like object
+        An open file to read from.
+
+    Returns
+    -------
+    grid_id : str
+        The ID of the Surfer ASCII grid.
+    shape : tuple = (n_northing, n_easting)
+        The number of grid points in the northing and easting dimension, respectively.
+    region : tuple = (west, east, south, north)
+        The grid region.
+    data_range : list = [min, max]
+        The minimum and maximum data values.
+
+    """
+    # DSAA is a Surfer ASCII GRD ID
+    grid_id = input_file.readline().strip()
+    # Read the grid shape (n_northing, n_easting)
+    shape = tuple(int(i.strip()) for i in input_file.readline().split())
+    # Our x points North, so the first thing we read north-south.
+    south, north = [float(i.strip()) for i in input_file.readline().split()]
+    west, east = [float(i.strip()) for i in input_file.readline().split()]
+    region = (west, east, south, north)
+    # The min and max of the data values (used for validation)
+    data_range = [float(i.strip()) for i in input_file.readline().split()]
+    return grid_id, shape, region, data_range
+
+
+def _check_surfer_integrity(field, shape, data_range):
+    """
+    Check that the grid matches the information from the header.
+    """
+    if field.shape != shape:
+        raise IOError(
+            "Grid shape {} doesn't match shape read from header {}.".format(
+                field.shape, shape
+            )
+        )
+    field_range = [field.min(), field.max()]
+    if not np.allclose(field_range, data_range):
+        raise IOError(
+            "Grid data range {} doesn't match range read from header {}.".format(
+                field_range, data_range
+            )
+        )


Test patch:
diff --git a/verde/tests/test_io.py b/verde/tests/test_io.py
new file mode 100644
index 0000000..046d17b
--- /dev/null
+++ b/verde/tests/test_io.py
@@ -0,0 +1,133 @@
+# pylint: disable=redefined-outer-name
+"""
+Test the I/O functions.
+"""
+import os
+from tempfile import NamedTemporaryFile
+from io import StringIO
+
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+from ..io import load_surfer, _read_surfer_header, _check_surfer_integrity
+
+
+@pytest.fixture
+def sample_grid():
+    "A simple grid for testing"
+    grid = StringIO(
+        """DSAA
+        2 6
+        1.0 1000.0
+        2.0 2000.0
+        11 61
+        11 21 31 41 51 61
+        12 22 32 42 52 27.5
+        """
+    )
+    return grid
+
+
+@pytest.fixture
+def sample_grid_nans():
+    "A simple grid with NaNs for testing"
+    grid = StringIO(
+        """DSAA
+        2 6
+        1.0 1000.0
+        2.0 2000.0
+        11 61
+        11 21 1.70141e38 41 1.70141e38 61
+        12 22 32 42 52 1.70141e38
+        """
+    )
+    return grid
+
+
+def test_load_surfer_header(sample_grid):
+    "Tests the surfer header parser"
+    grid_id, shape, region, data_range = _read_surfer_header(sample_grid)
+    assert grid_id == "DSAA"
+    assert shape == (2, 6)
+    npt.assert_allclose(region, (2, 2000, 1, 1000))
+    npt.assert_allclose(data_range, (11, 61))
+
+
+def test_load_surfer(sample_grid):
+    "Tests the main function on a StringIO file-like object"
+    grid = load_surfer(sample_grid)
+    assert grid.dims == ("northing", "easting")
+    assert grid.attrs["gridID"] == "DSAA"
+    assert grid.shape == (2, 6)
+    npt.assert_allclose(grid.northing.min(), 1)
+    npt.assert_allclose(grid.northing.max(), 1000)
+    npt.assert_allclose(grid.easting.min(), 2)
+    npt.assert_allclose(grid.easting.max(), 2000)
+    npt.assert_allclose(grid.min(), 11)
+    npt.assert_allclose(grid.max(), 61)
+    field = np.array([[11, 21, 31, 41, 51, 61], [12, 22, 32, 42, 52, 27.5]])
+    npt.assert_allclose(grid.values, field)
+
+
+def test_load_surfer_nans(sample_grid_nans):
+    "Tests the main function on a StringIO file-like object"
+    grid = load_surfer(sample_grid_nans)
+    assert grid.dims == ("northing", "easting")
+    assert grid.attrs["gridID"] == "DSAA"
+    assert grid.shape == (2, 6)
+    npt.assert_allclose(grid.northing.min(), 1)
+    npt.assert_allclose(grid.northing.max(), 1000)
+    npt.assert_allclose(grid.easting.min(), 2)
+    npt.assert_allclose(grid.easting.max(), 2000)
+    npt.assert_allclose(grid.min(), 11)
+    npt.assert_allclose(grid.max(), 61)
+    field = np.array([[11, 21, np.nan, 41, np.nan, 61], [12, 22, 32, 42, 52, np.nan]])
+    npt.assert_allclose(grid.values, field)
+
+
+def test_load_surfer_file(sample_grid):
+    "Tests the main function on a file on disk"
+    tmpfile = NamedTemporaryFile(delete=False)
+    tmpfile.close()
+    try:
+        with open(tmpfile.name, "w") as infile:
+            infile.write(sample_grid.getvalue())
+        grid = load_surfer(tmpfile.name)
+        assert grid.dims == ("northing", "easting")
+        assert grid.attrs["gridID"] == "DSAA"
+        assert grid.attrs["file"] == tmpfile.name
+        assert grid.shape == (2, 6)
+        npt.assert_allclose(grid.northing.min(), 1)
+        npt.assert_allclose(grid.northing.max(), 1000)
+        npt.assert_allclose(grid.easting.min(), 2)
+        npt.assert_allclose(grid.easting.max(), 2000)
+        npt.assert_allclose(grid.min(), 11)
+        npt.assert_allclose(grid.max(), 61)
+        field = np.array([[11, 21, 31, 41, 51, 61], [12, 22, 32, 42, 52, 27.5]])
+        npt.assert_allclose(grid.values, field)
+    finally:
+        os.remove(tmpfile.name)
+
+
+def test_load_surfer_integrity():
+    "Tests the surfer integrity check function"
+    field = np.array([[1, 2], [3, 4], [5, 6]])
+    _check_surfer_integrity(field, shape=(3, 2), data_range=(1, 6))
+    with pytest.raises(IOError):
+        _check_surfer_integrity(field, shape=(2, 3), data_range=(1, 6))
+    with pytest.raises(IOError):
+        _check_surfer_integrity(field, shape=(3, 2), data_range=(-2, 3))
+    # Make a corrupted grid file
+    wrong_shape = StringIO(
+        """DSAA
+        2 5
+        1.0 1000.0
+        2.0 2000.0
+        11 61
+        11 21 31 41 51 61
+        12 22 32 42 52 1.70141e38
+        """
+    )
+    with pytest.raises(IOError):
+        load_surfer(wrong_shape)