# 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)>
[54064 values with 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 ...
dy       float64 ...```

## 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()
```
```<matplotlib.collections.QuadMesh object at 0x7f19a4a806a0>
```

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()
```
```<matplotlib.collections.QuadMesh object at 0x7f199c9c04f0>
```

## 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.

```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 ...
dy       float64 ...
* 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 ...
dy       float64 ...
* 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()
```
```<matplotlib.collections.QuadMesh object at 0x7f199c869fc0>
```

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")
```
```[<matplotlib.lines.Line2D object at 0x7f199c86beb0>, <matplotlib.lines.Line2D object at 0x7f199c86bd90>, <matplotlib.lines.Line2D object at 0x7f199c86bdc0>, <matplotlib.lines.Line2D object at 0x7f199c73ecb0>, <matplotlib.lines.Line2D object at 0x7f199c73ec20>, <matplotlib.lines.Line2D object at 0x7f199c73ce50>, <matplotlib.lines.Line2D object at 0x7f199c73fcd0>, <matplotlib.lines.Line2D object at 0x7f199c73fdc0>, <matplotlib.lines.Line2D object at 0x7f199c73feb0>, <matplotlib.lines.Line2D object at 0x7f199c73ffa0>, <matplotlib.lines.Line2D object at 0x7f199c86be80>, <matplotlib.lines.Line2D object at 0x7f199c76c0d0>, <matplotlib.lines.Line2D object at 0x7f199c76c1c0>]
```

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()
```
```<matplotlib.collections.QuadMesh object at 0x7f199c66ca60>
```

## 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))
```
```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))
```
```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 ...
dy       float64 ...
* 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.694 seconds)

Gallery generated by Sphinx-Gallery