Source code for verde.io

# Copyright (c) 2017 The Verde Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Functions for input and output of grids in less common formats.
"""
import numpy as np
import xarray as xr


[docs] 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") # noqa: SIM115 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 ) )