Raster data and xarray#

Geospatial data primarily comes in two forms: raster data and vector data. This guide focuses on the first.

Raster data consists of rows and columns of rectangular cells. Their location in space is defined by the number of rows, the number of columns, a cell size along the rows, a cell size along the columns, the origin (x, y), and optionally rotation (x, y) – an affine matrix.

Typical examples of file formats containing raster data are:

  • GeoTIFF

  • ESRII ASCII

  • netCDF

  • IDF

In groundwater modeling, data commonly stored in raster format are:

  • Layer topology: the tops and bottoms of geohydrological layers

  • Layer properties: conductivity of aquifers and aquitards

  • Model output: heads or cell by cell flows

These data consist of values for every single cell. Xarray provides many conveniences for such data, such as plotting or selecting. To demonstrate, we’ll get some sample data provided with the imod package.

import xarray as xr
import imod

elevation = imod.data.ahn()["ahn"]
elevation
<xarray.DataArray 'ahn' (y: 218, x: 248)>
array([[ 5.195652e-01,  1.059200e+00,  6.309091e-01, ..., -1.719444e+00,
        -1.578696e+00, -1.501304e+00],
       [ 8.183333e-01,  7.568182e-01,  8.080000e-01, ..., -3.872800e+00,
        -1.928571e+00, -1.541667e+00],
       [ 8.673913e-01,  5.758824e-01, -1.428571e-03, ..., -5.064000e+00,
        -4.378750e+00, -3.507083e+00],
       ...,
       [-4.486000e+00, -5.040000e+00, -5.097600e+00, ..., -1.638400e+00,
        -1.656800e+00, -1.632400e+00],
       [-4.725600e+00, -4.526800e+00, -4.395652e+00, ..., -1.495600e+00,
        -1.671200e+00, -1.670000e+00],
       [-3.899600e+00, -4.262609e+00, -4.194400e+00, ..., -1.472400e+00,
        -1.599600e+00, -1.619600e+00]], dtype=float32)
Coordinates:
  * x        (x) float64 9.095e+04 9.105e+04 9.115e+04 ... 1.156e+05 1.156e+05
  * y        (y) float64 4.676e+05 4.674e+05 4.674e+05 ... 4.46e+05 4.458e+05
    dx       float64 100.0
    dy       float64 -100.0


Two dimensions: x, y#

This dataset represents some surface elevation in the west of the Netherlands, in the form of an xarray DataArray.

Xarray provides a “rich” representation of this data, note the x and y coordinates are shown above. We can use these coordinates for plotting, selecting, etc.

elevation.plot()
dx = 100.0, dy = -100.0

Out:

<matplotlib.collections.QuadMesh object at 0x7f40eab87f10>

This creates an informative plot.

We can also easily make a selection of a 10 by 10 km square, and plot the result:

selection = elevation.sel(x=slice(100_000.0, 110_000.0), y=slice(460_000.0, 450_000.0))
selection.plot()
dx = 100.0, dy = -100.0

Out:

<matplotlib.collections.QuadMesh object at 0x7f40e28a9120>

More dimensions#

Raster data can also be “stacked” to represent additional dimensions, such as height or time. Xarray is fully N-dimensional, and can directly represent these data.

Let’s start with a three dimensional example: a geohydrological layer model.

layermodel = imod.data.hondsrug_layermodel()
layermodel
<xarray.Dataset>
Dimensions:  (x: 500, y: 200, layer: 13)
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13
Data variables:
    top      (layer, y, x) float32 ...
    bottom   (layer, y, x) float32 ...
    k        (layer, y, x) float32 ...
    idomain  (layer, y, x) float32 ...


This dataset contains multiple variables. We can take a closer look at the the “top” variable, which represents the top of every layer.

top = layermodel["top"]
top
<xarray.DataArray 'top' (layer: 13, y: 200, x: 500)>
[1300000 values with dtype=float32]
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13


This DataArray has three dimensions: layer, y, x. We can’t make a planview plot of this entire dataset: every (x, y) locations has as many values as layers. A single layer can be selected and shown as follows:

top.sel(layer=1).plot()
dx = 25.0, dy = -25.0, layer = 1

Out:

<matplotlib.collections.QuadMesh object at 0x7f40e2961900>

Xarray doesn’t favor specific dimensions. We can select a value along the y-axis just as easily, to create a cross-section.

section = top.sel(y=560_000.0, method="nearest")
section.plot.line(x="x")
y = 5.6e+05, dx = 25.0, dy = -25.0

Out:

[<matplotlib.lines.Line2D object at 0x7f40e28406d0>, <matplotlib.lines.Line2D object at 0x7f40e2840730>, <matplotlib.lines.Line2D object at 0x7f40e2840850>, <matplotlib.lines.Line2D object at 0x7f40e2840970>, <matplotlib.lines.Line2D object at 0x7f40e2840a90>, <matplotlib.lines.Line2D object at 0x7f40e2840bb0>, <matplotlib.lines.Line2D object at 0x7f40e2840cd0>, <matplotlib.lines.Line2D object at 0x7f40e2840df0>, <matplotlib.lines.Line2D object at 0x7f40e2840f10>, <matplotlib.lines.Line2D object at 0x7f40e2841030>, <matplotlib.lines.Line2D object at 0x7f40e2841150>, <matplotlib.lines.Line2D object at 0x7f40e2840700>, <matplotlib.lines.Line2D object at 0x7f40e2841360>]

Xarray provides us a lot of convenience compared to working with traditional two dimensional rasters: rather than continuously loop over the data of single timesteps or layers, we can process them in a single command.

For example, to compute the thickness of every layer:

thickness = layermodel["top"] - layermodel["bottom"]
thickness
<xarray.DataArray (layer: 13, y: 200, x: 500)>
array([[[ 0.        ,  0.        ,  0.        , ...,  3.8964    ,
          3.8964    ,  3.8964    ],
        [ 0.        ,  0.        ,  0.        , ...,  3.8964    ,
          3.8964    ,  3.8964    ],
        [ 0.        ,  0.        ,  0.        , ...,  3.8964    ,
          3.8964    ,  3.8964    ],
        ...,
        [ 1.302     ,  1.302     ,  1.302     , ...,  0.        ,
          0.        ,  0.        ],
        [ 1.302     ,  1.302     ,  1.302     , ...,  0.        ,
          0.        ,  0.        ],
        [ 1.302     ,  1.302     ,  1.302     , ...,  0.        ,
          0.        ,  0.        ]],

       [[ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
...
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ]],

       [[63.152397  , 63.152397  , 63.152397  , ..., 94.64      ,
         94.64      , 94.64      ],
        [63.152397  , 63.152397  , 63.152397  , ..., 94.64      ,
         94.64      , 94.64      ],
        [63.152397  , 63.152397  , 63.152397  , ..., 94.64      ,
         94.64      , 94.64      ],
        ...,
        [74.6828    , 74.6828    , 74.6828    , ..., 78.1022    ,
         78.1022    , 78.1022    ],
        [74.6828    , 74.6828    , 74.6828    , ..., 78.1022    ,
         78.1022    , 78.1022    ],
        [74.6828    , 74.6828    , 74.6828    , ..., 78.1022    ,
         78.1022    , 78.1022    ]]], dtype=float32)
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13


This is easily multiplied, then summed over the layer dimensions to provide us a map of the total transmissivity:

transmissivity = layermodel["k"] * thickness
total = transmissivity.sum("layer")
total.plot()
dx = 25.0, dy = -25.0

Out:

<matplotlib.collections.QuadMesh object at 0x7f40e26f3a00>

Input and output#

The imod package started as a collection of functions to read IDF files into xarray DataArrays. By convention, IDF files store the coordinates of the extra dimensions (layer, time) in the file name. imod.idf.save() will automatically generate these names from a DataArray.

Let’s demonstrate by writing the transmissivity computed above to IDF. (We’ll do this in a temporary directory to keep things tidy.)

tempdir = imod.util.temporary_directory()
imod.idf.save(tempdir / "transmissivity", transmissivity)

Note

tempdir is a Python pathlib.Path object. These objects are very convenient for working with paths; we can easily check if paths exists, join paths with /, etc.

Let’s check which files have been written in the temporary directory:

filenames = [path.name for path in tempdir.iterdir()]
print("\n".join(filenames))

Out:

transmissivity_l10.idf
transmissivity_l12.idf
transmissivity_l13.idf
transmissivity_l6.idf
transmissivity_l9.idf
transmissivity_l4.idf
transmissivity_l5.idf
transmissivity_l3.idf
transmissivity_l11.idf
transmissivity_l1.idf
transmissivity_l2.idf
transmissivity_l7.idf
transmissivity_l8.idf

Just as easily, we can read all IDFs back into a single DataArray. We can do so by using a wildcard. These wildcards function according to the rules of Glob via the python glob module.

Note that every IDF has to have identical x-y coordinates: files with different cell sizes or extents will not be combined automatically.

back = imod.idf.open(tempdir / "*.idf")
back
<xarray.DataArray 'transmissivity' (layer: 13, y: 200, x: 500)>
dask.array<stack, shape=(13, 200, 500), dtype=float32, chunksize=(1, 200, 500), chunktype=numpy.ndarray>
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int64 1 2 3 4 5 6 7 8 9 10 11 12 13


These glob patterns are quite versatile, and may be used to filter as well.

selection = imod.idf.open(tempdir / "transmissivity_l[1-5].idf")
selection
<xarray.DataArray 'transmissivity' (layer: 5, y: 200, x: 500)>
dask.array<stack, shape=(5, 200, 500), dtype=float32, chunksize=(1, 200, 500), chunktype=numpy.ndarray>
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int64 1 2 3 4 5


Rather commonly, the paths of the IDFs are not named according to consistent rules. In such cases, we can manually specify how the name should be interpreted via the pattern argument.

back = imod.idf.open(tempdir / "*.idf", pattern="{name}_l{layer}")
back
<xarray.DataArray 'transmissivity' (layer: 13, y: 200, x: 500)>
dask.array<stack, shape=(13, 200, 500), dtype=float32, chunksize=(1, 200, 500), chunktype=numpy.ndarray>
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int64 1 2 3 4 5 6 7 8 9 10 11 12 13


See the documenation of imod.idf.open() for more details.

Other raster formats#

IDF is one raster format, but there are many more. imod.rasterio.open() wraps the rasterio Python package (which in turn wraps GDAL) to provide access to many GIS raster formats. imod.rasterio.open() and imod.rasterio.save() work exactly the same as the respective IDF functions, except they write to a different format.

For example, to write the transsmisivity to GeoTIFF:

imod.rasterio.save(tempdir / "transmissivity.tif", transmissivity)
filenames = [path.name for path in tempdir.iterdir()]
print("\n".join(filenames))

Out:

transmissivity_l4.tif
transmissivity_l10.idf
transmissivity_l12.idf
transmissivity_l13.idf
transmissivity_l6.idf
transmissivity_l9.idf
transmissivity_l11.tif
transmissivity_l8.tif
transmissivity_l4.idf
transmissivity_l5.idf
transmissivity_l5.tif
transmissivity_l6.tif
transmissivity_l12.tif
transmissivity_l3.idf
transmissivity_l9.tif
transmissivity_l7.tif
transmissivity_l11.idf
transmissivity_l13.tif
transmissivity_l1.idf
transmissivity_l3.tif
transmissivity_l2.tif
transmissivity_l1.tif
transmissivity_l10.tif
transmissivity_l2.idf
transmissivity_l7.idf
transmissivity_l8.idf

Note imod.rasterio.save() will split the extension off the path, infer the GDAL driver, attach the additional coordinates to the file name, and re-attach the extension.

netCDF#

The final format to discuss here is netCDF. Two dimensional raster files are convenient for viewing, as every file corresponds with a single “map view” in a GIS viewer. However, they are not convenient for storing many timesteps or many layers: especially long running simulations will generate hundreds, thousands, or even millions of files.

netCDF is a file format specifically designed for multi-dimensional data, and allows us to conveniently bundle all data in a single file. Xarray objects can directly be written to netCDF, as the data model of xarray itself is based on the netCDF data model.

With netCDF, there is no need to encode the different dimension labels in the the name: they are stored directly in the file instead.

layermodel.to_netcdf(tempdir / "layermodel.nc")
back = xr.open_dataset(tempdir / "layermodel.nc")
back
<xarray.Dataset>
Dimensions:  (x: 500, y: 200, layer: 13)
Coordinates:
  * x        (x) float64 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
  * y        (y) float64 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
    dx       float64 25.0
    dy       float64 -25.0
  * layer    (layer) int32 1 2 3 4 5 6 7 8 9 10 11 12 13
Data variables:
    top      (layer, y, x) float32 ...
    bottom   (layer, y, x) float32 ...
    k        (layer, y, x) float32 ...
    idomain  (layer, y, x) float32 ...


Coordinate reference systems (CRS)#

Reprojection from one CRS to another is a common frustration. Since the data in an xarray DataArray is always accompanied by its x and y coordinates, we can easily reproject the data. See the examples.

Total running time of the script: ( 0 minutes 1.782 seconds)

Gallery generated by Sphinx-Gallery