Input/Output (IO)

This section documents various methods you can use to bring in data and write output to a file.

• Input and output of Landlab raster-grid data in ESRI Ascii Raster format
  • ESRI ASCII functions
  • BadHeaderLineError
  • DataSizeError
  • Error
  • KeyTypeError
  • KeyValueError
  • MismatchGridDataSizeError
  • MismatchGridXYLowerLeft
  • MismatchGridXYSpacing
  • MissingRequiredKeyError
  • read_asc_header
  • read_esri_ascii
  • write_esri_ascii
• Input and output in native landlab format
  • Read Landlab native
  • load_grid
  • save_grid
• Read and Write netCDF Files
  • Reading
  • Writing
  • Error Handling
• Shapefiles
  • Reading
• Output of Landlab raster- and hex-grid data in Wavefront OBJ format
  • OBJ functions
  • write_obj
• Output of Landlab grid and field data in Legacy VTK format
  • write_legacy_vtk

Module contents

Modules that read/write ModelGrids from various file formats.

exception BadHeaderLineError(line)

Bases: Error

Raise this error for a bad header is line.

exception DataSizeError(size, expected_size)

Bases: Error

Raise this error if the size of data does not match the header.

exception KeyTypeError(key, expected_type)

Bases: Error

Raise this error when a header’s key value is of the wrong type.

exception KeyValueError(key, message)

Bases: Error

Raise this error when a header’s key value has a bad value.

exception MismatchGridDataSizeError(size, expected_size)

Bases: Error

Raise this error if the data size does not match the grid size.

exception MismatchGridXYLowerLeft(llc, expected_llc)

Bases: Error

Raise this error if the file lower left does not match the grid.

exception MismatchGridXYSpacing(dx, expected_dx)

Bases: Error

Raise this error if the file cell size does not match the grid dx.

exception MissingRequiredKeyError(key)

Bases: Error

Raise this error when a header is missing a required key.

read_asc_header(asc_file)

Read header information from an ESRI ASCII raster file.

The header contains the following variables,

• ncols: Number of cell columns

• nrows: Number of cell rows

• xllcenter or xllcorner: X (column) coordinate of lower-left

  coordinate of grid (by center or lower-left corner of the cell)

• yllcenter, yllcorner: Y (row) coordinate of lower-left

  coordinate of grid (by center or lower-left corner of the cell)

• cellsize: Grid spacing between rows and columns

• nodata_value: No-data value (optional)

Parameters:

asc_file (file_like) – File-like object from which to read header.

Returns:

Header as key-value pairs.

Return type:

dict

Raises:

• MissingRequiredKeyError – The header is missing a required key.

• KeyTypeError – The header has the key but its values is of the wrong type.

Examples

>>> from io import StringIO
>>> from landlab.io.esri_ascii import read_asc_header

>>> contents = '''
... nrows 100
... ncols 200
... cellsize 1.5
... xllcenter 0.5
... yllcenter -0.5
... '''

>>> hdr = read_asc_header(StringIO(contents))
>>> hdr["nrows"], hdr["ncols"]
(100, 200)
>>> hdr["cellsize"]
1.5
>>> hdr["xllcenter"], hdr["yllcenter"]
(0.5, -0.5)

MissingRequiredKeyError is raised if the header does not contain all of the necessary keys.

>>> contents = '''
... ncols 200
... cellsize 1.5
... xllcenter 0.5
... yllcenter -0.5
... '''
>>> read_asc_header(StringIO(contents)) 
Traceback (most recent call last):
MissingRequiredKeyError: nrows

KeyTypeError is raised if a value is of the wrong type. For instance, nrows and ncols must be int.

>>> contents = '''
... nrows 100.5
... ncols 200
... cellsize 1.5
... xllcenter 0.5
... yllcenter -0.5
... '''
>>> read_asc_header(StringIO(contents)) 
Traceback (most recent call last):
KeyTypeError: Unable to convert nrows to <type 'int'>

read_esri_ascii(asc_file, grid=None, reshape=False, name=None, halo=0)

Read RasterModelGrid from an ESRI ASCII file.

Read data from asc_file, an ESRI ASCII file, into a RasterModelGrid. asc_file is either the name of the data file or is a file-like object.

The grid and data read from the file are returned as a tuple (grid, data) where grid is an instance of RasterModelGrid and data is a numpy array of doubles with that has been reshaped to have the number of rows and columns given in the header.

Parameters:

• asc_file (str of file-like) – Data file to read.

• reshape (boolean, optional) – Reshape the returned array, otherwise return a flattened array.

• name (str, optional) – Add data to the grid as a named field.

• grid (grid , optional) – Adds data to an existing grid instead of creating a new one.

• halo (integer, optional) – Adds outer border of depth halo to the grid.

Returns:

(grid, data) – A newly-created RasterModelGrid and the associated node data.

Return type:

tuple

Raises:

• DataSizeError – Data are not the same size as indicated by the header file.

• MismatchGridDataSizeError – If a grid is passed, and the size of the grid does not agree with the size of the data.

• MismatchGridXYSpacing – If a grid is passed, and the cellsize listed in the heading does not match the node spacing of the grid.

• MismatchGridXYLowerLeft – If a grid is passed and the xllcorner and yllcorner do not match that of the grid.

Examples

>>> from landlab.io import read_esri_ascii
>>> from io import StringIO

>>> contents = '''
... ncols         3
... nrows         4
... xllcorner     1.
... yllcorner     2.
... cellsize      10.
... NODATA_value  -1
... 0. 1. 2.
... 3. 4. 5.
... 6. 7. 8.
... 9. 10. 11.
... '''

>>> (grid, data) = read_esri_ascii(StringIO(contents))

The returned grid is a RasterModelGrid with 4 rows and 3 columns.

>>> grid
RasterModelGrid((4, 3), xy_spacing=(10.0, 10.0), xy_of_lower_left=(1.0, 2.0))

Note that the first row of values is the bottom-most of the data file.

>>> data.reshape(grid.shape)
array([[  9.,  10.,  11.],
       [  6.,   7.,   8.],
       [  3.,   4.,   5.],
       [  0.,   1.,   2.]])

>>> (grid, data) = read_esri_ascii(StringIO(contents), halo=1)

Because of the halo, the returned grid now has two more rows and columns than before.

>>> grid
RasterModelGrid((6, 5), xy_spacing=(10.0, 10.0), xy_of_lower_left=(-9.0, -8.0))
>>> data.reshape(grid.shape)
array([[ -1.,  -1.,  -1.,  -1.,  -1.],
       [ -1.,   9.,  10.,  11.,  -1.],
       [ -1.,   6.,   7.,   8.,  -1.],
       [ -1.,   3.,   4.,   5.,  -1.],
       [ -1.,   0.,   1.,   2.,  -1.],
       [ -1.,  -1.,  -1.,  -1.,  -1.]])

read_shapefile(file, dbf=None, store_polyline_vertices=True, points_shapefile=None, points_dbf=None, link_fields=None, node_fields=None, link_field_conversion=None, node_field_conversion=None, link_field_dtype=None, node_field_dtype=None, threshold=0.0)

Read shapefile and create a NetworkModelGrid.

There are a number of assumptions that are required about the shapefile.

• The shapefile must be a polyline shapefile.

• All polylines must be their own object (e.g. no multi-part polylines).

• Polyline endpoints match perfectly.

You might notice that there is no write_shapefile function. If this is something you need for your work, please make a GitHub issue to start this process.

Parameters:

• file (str or file-like) – File path or file-like of a valid polyline shapefile

• dbf (file-like, optional) – If file is file-like, the dbf must also be passed.

• store_polyline_vertices (bool, optional) – If True (default), store the vertices of the polylines in the at_link fields x_of_polyline and y_of_polyline.

• points_shapefile (str or file-like) – File path or file-like of a valid point shapefile.

• points_dbf (file-like, optional) – If file is file-like, the dbf must also be passed.

• link_fields (list, optional) – List of polyline shapefile attributes to import as landlab at-link fields. Default is to import all.

• node_fields (list, optional) – List of point shapefile attributes to import as landlab at-node fields. Default is to import all.

• link_field_conversion (dict, optional) – Dictionary mapping polyline shapefile field names to desired at link field names. Default is no remapping.

• node_field_conversion (dict, optional) – Dictionary mapping node shapefile field names to desired at node field names. Default is no remapping.

• link_field_dtype (dict, optional) – Dictionary mapping node shapefile field names to desired dtype. Default is no change to dtype.

• node_field_dtype (dict, optional) – Dictionary mapping node shapefile field names to desired dtype. Default is no change to dtype.

• threshold (float, optional) – Maximum distance between a point in the point shapefile and a polyline junction in the polyline shapefile. Units are the same as in the shapefiles. Default is zero (requiring perfect overlap).

Returns:

grid – The network model grid will have nodes at the endpoints of the polylines, and links that connect these nodes. Any fields associated with the shapefile will be added as at-link fields. If a point shapefile is provided those values will be added as at-node fields.

Return type:

NetworkModelGrid instance

Examples

First, we make a simple shapefile

>>> from io import BytesIO
>>> import os

>>> import shapefile
>>> shp = BytesIO()
>>> shx = BytesIO()
>>> dbf = BytesIO()
>>> w = shapefile.Writer(shp=shp, shx=shx, dbf=dbf)
>>> w.shapeType = shapefile.POLYLINE
>>> w.field("spam", "N")
>>> w.line([[[5,5],[10,10]]])
>>> w.record(37)
>>> w.line([[[5,0],[5,5]]])
>>> w.record(100)
>>> w.line([[[5,5],[0,10]]])
>>> w.record(239)
>>> w.close()

Now create a NetworkModelGrid with read_shapefile:

>>> from landlab.io import read_shapefile
>>> grid = read_shapefile(shp, dbf=dbf)
>>> grid.nodes
array([0, 1, 2, 3])
>>> grid.x_of_node
array([  5.,   5.,   0.,  10.])
>>> grid.y_of_node
array([  0.,   5.,  10.,  10.])
>>> grid.nodes_at_link
array([[0, 1],
       [2, 1],
       [1, 3]])
>>> assert "spam" in grid.at_link
>>> grid.at_link["spam"]
array([100, 239, 37])

Next lets also include a points file. First create both shapefiles.

>>> shp = BytesIO()
>>> shx = BytesIO()
>>> dbf = BytesIO()
>>> w = shapefile.Writer(shp=shp, shx=shx, dbf=dbf)
>>> w.shapeType = shapefile.POLYLINE
>>> w.field("spam", "N")
>>> w.line([[[5,5],[10,10]]])
>>> w.record(37)
>>> w.line([[[5,0],[5,5]]])
>>> w.record(100)
>>> w.line([[[5,5],[0,10]]])
>>> w.record(239)
>>> w.close()

>>> p_shp = BytesIO()
>>> p_shx = BytesIO()
>>> p_dbf = BytesIO()
>>> p_w = shapefile.Writer(shp=p_shp, shx=p_shx, dbf=p_dbf)
>>> p_w.shapeType = shapefile.POINT
>>> p_w.field("eggs", "N")
>>> p_w.point(5, 0)
>>> p_w.record(2)
>>> p_w.point(5, 5)
>>> p_w.record(4)
>>> p_w.point(0, 10)
>>> p_w.record(8)
>>> p_w.point(10, 10)
>>> p_w.record(6)
>>> p_w.close()

Now read in both files together.

>>> grid = read_shapefile(shp, dbf=dbf, points_shapefile=p_shp, points_dbf=p_dbf)
>>> grid.nodes
array([0, 1, 2, 3])
>>> grid.x_of_node
array([  5.,   5.,   0.,  10.])
>>> grid.y_of_node
array([  0.,   5.,  10.,  10.])
>>> grid.nodes_at_link
array([[0, 1],
       [2, 1],
       [1, 3]])
>>> assert "spam" in grid.at_link
>>> grid.at_link["spam"]
array([100, 239, 37])
>>> assert "eggs" in grid.at_node
>>> grid.at_node["eggs"]
array([2, 4, 8, 6])

write_esri_ascii(path, fields, names=None, clobber=False)

Write landlab fields to ESRI ASCII.

Write the data and grid information for fields to path in the ESRI ASCII format.

Parameters:

• path (str) – Path to output file.

• fields (field-like) – Landlab field object that holds a grid and associated values.

• names (iterable of str, optional) – Names of the fields to include in the output file. If not provided, write all fields.

• clobber (boolean) – If path exists, clobber the existing file, otherwise raise an exception.

Examples

>>> import numpy as np
>>> import os
>>> from landlab import RasterModelGrid
>>> from landlab.io.esri_ascii import write_esri_ascii

>>> grid = RasterModelGrid((4, 5), xy_spacing=(2., 2.))
>>> grid.at_node["air__temperature"] = np.arange(20.0)
>>> files = write_esri_ascii("test.asc", grid)  
>>> [os.path.basename(name) for name in sorted(files)])  
['test.asc']

>>> _ = grid.add_field("land_surface__elevation", np.arange(20.), at="node")
>>> grid.at_node["land_surface__elevation"] = np.arange(20.0)
>>> files = write_esri_ascii("test.asc", grid))  
>>> [os.path.basename(name) for name in sorted(files)])  
['test_air__temperature.asc', 'test_land_surface__elevation.asc']

write_legacy_vtk(path, grid, z_at_node='topographic__elevation', fields=None, clobber=False)

Write grid and field to a legacy VTK format file or file-like object.

Parameters:

• path (file-like) – Path to file or a file-like object

• grid (Landlab grid object) – The grid for which to output data

• z_at_node (str or (n_nodes, ) ndarray) – Field name or array of values to use for z coordinate

• fields (list of str (optional)) – List of node fields to output; default is all node fields

• clobber (bool (optional)) – Ok to overwrite existing file (default False)

Examples

>>> import io
>>> import numpy as np
>>> from landlab import HexModelGrid
>>> from landlab.io.legacy_vtk import write_legacy_vtk

>>> grid = HexModelGrid((3, 2))
>>> topo = grid.add_zeros("topographic__elevation", at="node")
>>> topo[:] = np.arange(len(topo))
>>> water = grid.add_zeros("surface_water__depth", at="node")
>>> water[:] = 0.1 * (7.0 - topo)

>>> vtk_file = write_legacy_vtk(io.StringIO(), grid)
>>> lines = vtk_file.getvalue().splitlines()
>>> print(lines[0])
# vtk DataFile Version 2.0
>>> for i in range(5, 13):
...     print(lines[i])
POINTS 7 float
0.5 0.0 0.0
1.5 0.0 1.0
0.0 0.866025 2.0
1.0 0.866025 3.0
2.0 0.866025 4.0
0.5 1.732051 5.0
1.5 1.732051 6.0
>>> for i in range(14, 21):
...     print(lines[i])
CELLS 6 24
3 3 0 1
3 3 2 0
3 4 3 1
3 5 2 3
3 6 3 4
3 6 5 3
>>> for i in range(22, 29):
...     print(lines[i])
CELL_TYPES 6
5
5
5
5
5
5
>>> for i in range(30, 49):
...     print(lines[i])
POINT_DATA 7
SCALARS topographic__elevation float 1
LOOKUP_TABLE default
0.0
1.0
2.0
3.0
4.0
5.0
6.0
SCALARS surface_water__depth float 1
LOOKUP_TABLE default
0.7
0.6
0.5
0.4
0.3
0.2
0.1

write_obj(path, grid, field_for_z='topographic__elevation', clobber=False)

Write landlab grid and one field to Wavefront OBJ.

Parameters:

• path (str, or file-like) – Path to output file.

• grid (Landlab grid object) – Landlab grid object that includes associated values.

• field_for_z (str, optional) – Name of a field to use for the z-values of the OBJ file.

• clobber (boolean, optional) – If path exists, clobber the existing file, otherwise raise an exception.

Returns:

The input path used to write the OBJ file.

Return type:

str or file-like

Examples

>>> import io
>>> from landlab import HexModelGrid, RasterModelGrid
>>> from landlab.io.obj import write_obj

>>> grid = HexModelGrid((3, 2), spacing=2.)
>>> z = grid.add_zeros("topographic__elevation", at="node")
>>> z[3] = 1.0

>>> obj_file = write_obj(io.StringIO(), grid)
>>> print(obj_file.getvalue())
# landlabgrid
#
g landlabgrid
v 1.0 0.0 0.0
v 3.0 0.0 0.0
v 0.0 1.732051 0.0
v 2.0 1.732051 1.0
v 4.0 1.732051 0.0
v 1.0 3.464102 0.0
v 3.0 3.464102 0.0
f 4// 1// 2//
f 4// 3// 1//
f 5// 4// 2//
f 6// 3// 4//
f 7// 4// 5//
f 7// 6// 4//