Raster (Grid) Handling¶

This section introduces geospatial analysis of raster (gridded) data with gdal and rasterstats. Download sample raster datasets from River Architect. In particular, this section uses GeoTIFF rasters located in RiverArchitect/SampleData/01_Conditions/2100_sample/.

Requirements:

• Make sure to understand gridded raster data, in particular, the .tif (GeoTIFF) format
• Accomplish the QGIS tutorial

Tips:

• While OSGeo's ogr module is useful for shapefile handling, raster data are best handled by gdal.
• The functions featured in this section are introduced are partially also implemented in flusstools.
• To use those functions, make sure flusstools is installed and import it as follows: from flusstools import geotools. Some of the functions shown in this tutorial can then be used with geotools.function_name().

Load raster¶

Open Existing Raster Data¶

Raster data can be opened in Python code as an instance of gdal.Open("FILENAME"). The following code block provides a function to open any raster specified with the file_name input argument. One of the most important elements when dealing with raster data is the raster band, which takes on a similar data carrier role as GetLayer in shapefile handling. To access the raster band, the below-shown open_raster function:

1. Enables error and warning feedback with gdal.UseExceptions() (this step is vital when using gdal).
2. Opens the provided raster file_name embraced by try - except statements to inform if and why a potential error occurred while opening the raster.
3. Opens the raster band number stated in the optional band_number keyword argument with raster_band = raster.GetRasterBand(band_number) (the default value is 1).
4. Returns the raster and raster band objects.

Tip: This function is available in flusstools. For instance, use it with from flusstools import geotools and geotools.open_raster("path/to/raster.tif").

To use the open_raster function, call it with a file name as shown in the following code block with the h001000.tif raster from the River Architect sample data. The script immediately closes the raster again by overwriting the variable with a None-instance to avoid the GeoTIFF file being locked afterward.

Raster Band Statistics and Toolbox Scripts¶

Once the raster and its band are loaded with the above open_raster function, we can access statistical information (e.g., the minimum or the maximum), identify the no-data value (i.e., a pre-defined value that is assigned to pixels without value), or the type of units used.

Python scripts for processing geospatial data can also be embedded as plugins in GIS desktop applications (e.g., as plugins in QGIS or Toolbox in ArcGIS Pro). To run a Python script in a GIS desktop application, it should be written as a standalone script that can receive input arguments. The interested reader can learn more about implementing plugins in QGIS in the QGIS docs.

Note: QGIS wraps many external algorithm types, which are available through the QGIS Processing Toolbox. These algorithms belong, for example, to SAGA or GRASS GIS.

Here we will only write the next code block so that it can be run in a console/terminal application as a standalone script (recall the instructions for writing standalone script).

To run this script, save it as raster_band_info.py (e.g., in C:\temp in Windows or ~/temp/ in Linux) and navigate to the script directory in a terminal interface (e.g., in PyCharm's Terminal, Anaconda Prompt, or the Linux Terminal) using the cd command. Then run the script to get statistics of the water depth raster h001000.tif with (in Windows):

C:\temp\ python raster_band_info.py 1 "C:\temp\geodata\rasters\h001000.tif"

On Linux and with the vflussenv activated tap, for example:

python raster_band_info.py 1 "~/temp/geodata/rasters/h001000.tif"

Band minimum:  0.0
Band maximum:  7.0613012313843
No-data value:  -3.4028234663852886e+38
Band unit type:
Band has no ColorTable.
ColorTable: None

Create and Save a Raster (from Array)¶

Raster Drivers¶

Just like for shapefile files, the appropriate gdal driver (analogous to ogr drivers) must be loaded to save a raster. To get a full list of gdal raster drivers run:

Raster Data Types¶

The output raster pixels can be one of the following data types (source: gdal.org/doxygen/):

• GDT_Unknown Unknown or unspecified type
• GDT_Byte 8 bit unsigned integer
• GDT_UInt16 16 bit unsigned integer
• GDT_Int16 16 bit signed integer
• GDT_UInt32 32 bit unsigned integer
• GDT_Int32 32 bit signed integer
• GDT_Float32 32 bit floating point
• GDT_Float64 64 bit floating point
• GDT_CInt16 Complex Int16
• GDT_CInt32 Complex Int32
• GDT_CFloat32 Complex Float32
• GDT_CFloat64 Complex Float64

Create a Raster (Array to Raster)¶

Knowing the basics of raster handling, data types, and Python, we can create a raster from a numeric array. Since a raster is basically a georeferenced array, it is convenient to convert a numpy array into a raster (band). The function blocks feature the conversion of a numpy array into a GeoTIFF raster according to the following workflow:

1. Check out the GeoTIFF driver (driver = gdal.GetDriverByName('GTiff')).
2. Retrieve the array size and (number of rows rows and columns cols).
3. Create a new GeoTIFF raster (new_raster = driver.Create(file_name, cols, rows, 1, eType=rdtype)), where
   • file_name is the directory and name of the new raster file ending on .tif (e.g., "C:\\temp\\rasters\\new.tif").
   • cols, rows represent the array shape, and eType is the geospatial data type (see above list)
4. Set the geographic origin, which is stored in the origin (tuple) parameter and defines the pixel_width and pixel_height (pixel units defined with srs - see below).
5. Replace np.nan values in the numpy array with nan_value.
6. Instantiate a band object, set the NoDataValue to nan_value, and write the array to the band.
7. Create a spatial reference system object (srs) as a function of the epsg input parameter and export it to WKT format.
8. Release the raster (flush from cache).

Note: The units defined with the epsg number drive the pixel size, where pixel_width and pixel_height are multipliers of that unit. In the case of epsg=3857, the unit is meters and pixel_width=10 combined with pixel_height=20 creates 10-m wide and 20-m high pixels. In the case of epsg=4326, the unit is (geographic) degrees and 1 degree by 1-degree pixels can have the size of a count(r)y.

To call the function for writing a random numpy array, we can now use the create_raster() function (also available from geotools.create_raster()):

Raster Calculus (Raster / Band to Array)¶

The procedure described with the create_raster() function can be used vice-versa, to create numpy array from raster bands. The raster converted to a numpy array enables algebraic or other logical operations to be applied to existing raster data.

Need an example? In the RiverArchitect SampleData, the units of the water depth raster h001000.tif are in U.S. customary feet and the units of the flow velocity raster u001000.tif are in feet per second. However, to calculate the Froude number (involves the gravity constant) for every pixel based on the two rasters (water depth and flow velocity), it is convenient to convert both rasters into m and m/s, respectively. For this purpose, the following code block features another re-usable, custom function that loads a raster as an array and overwrites NoDataValues with np.nan (raster and band can be instantiated with the above open_raster function):

The raster2array() function is also included in flusstools: geotools.raster2array()

Challenge: The raster2array function returns a tuple, where output[0] corresponds to the array and output[1] is the geo-transformation. Can you optimize the way how this information is returned?

Finally, to create a Froude number GeoTIFF raster, the following code block makes use of the raster2array function for converting the water depth and flow velocity GeoTIFF rasters into a numpy array, performing simple algebraic calculations to convert the rasters to m and m/s, respectively, and save the resulting GeoTIFF file. In detail, the workflow involves to:

• Define the input raster file names with directories (h_file and u_file),
• Load original rasters as ndarray with the raster2array() function and get the original GeoTransform description,
• Convert all values from U.S. customary feet to S.I. metric units (recall the feet_to_meter() function from the Python basics), and
• Save a new copy of the raster.

Reproject a Raster¶

The transformation (and reprojection) of a raster into another coordinate system involves rotating, shifting, and shearing pixels. If one of these operations is skipped, the reprojected raster may be squeezed, twisted, or placed anywhere in the world but not where it should be placed. Therefore, the approach for the reprojection of a raster into another coordinate reference system implies the following steps:

1. Retrieve the source and target spatial reference systems (e.g., derive from a gdal.Dataset or an EPSG authority code).
2. Read the geo transformation of the source dataset ( gdal.Dataset.GetGeoTransform()).
3. Derive the number of pixels and the spacing between pixels in the new (reprojected) dataset.
4. Instantiate the new (reprojected) dataset.
5. Project an image of the source dataset onto the new (reprojected) dataset (gdal.ReprojectImage()).

The spatial reference system can be derived from a dataset with the explanations in the shapefile section by writing a get_srs() function. The following code block shows the get_srs() function (uses the osr library from osgeo / gdal ), which is also integrated in flusstools (geotools.get_srs()).

With the open_raster() and get_srs() functions, we have all necessary ingredients to accomplish the raster reprojection workflow in another funcation called reproject_raster(). An additional feature of the function is that it ensures the correct use of osr.CoordinateTransformation, which behaves differently under gdal 3.0 compared with older gdal versions (read on OSGeo's GitHub page).

Using the reproject_raster function in a Python script requires a source dataset and another (orientation) dataset with the new coordinate system into which the source dataset will be projected. The following example shows how to re-project the above-created Froude-number raster into the EPSG=3857 coordinate system for viewing it in QGIS on a Google Satellite basemap (recall basemaps in QGIS tutorial). As orientation data set, use web_frame.tif, which was created on the Google Satellite basemap.

With the get_srs() function that automatically detects the raster projection and spatial reference system, we can use the create_raster() function to reproject the above-created Fr1000cfs.tif raster (e.g., to epsg=4326).

Plotted in QGIS, the reprojected Froude-number raster looks like this:

The reproject_raster function is also available (slightly modified) in flusstools (geotools.reproject_raster()), where saving the new reprojected raster is embedded in the function (automatically appends the syllable "_epsg[NO]" to the original file name).

Note: To display multiple rasters with different coordinate systems on the same map in QGIS, the coordinate systems must be harmonized in most cases. While QGIS will automatically propose to convert a raster with a non-project-compliant CRS to the project CRS, it also has a dedicated function for adjusting raster coordinate systems: In QGIS, click on the Raster menu > Projections > Warp (Reproject).... Select the raster(s) to reproject (i.e., the raster(s) to harmonize with the project coordinate system). However, Warp may not perform all reprojection steps as desired and lead to wrong placements of the new raster. The Warp method is also available in Python through gdal.Warp (read the docs):

kwargs={'format': 'GTiff', 'geoloc': True}
gdal.Warp(TARGET_GEO_TIFF_FILE_NAME, SOURCE_GEO_TIFF_FILE_NAME, **kwargs)

The coordinate transformation fails when no transformation between the indicated source and target spatial reference system can be established (i.e., gdal does not know the transformation). This problem occurs often when old, regional coordinate systems are transformed into coordinate systems for web applications (e.g., EPSG=3857). Read more in the gdal docs.`

Zonal Statistics for Morphological Units¶

In hydraulic and geospatial analyses, the question of statistical values of certain areas of one or more rasters often arises. For instance, we may be interested in mean values and standard deviations in specific zones of a water body. For this purpose, Zonal statistics enable the delineation of an area of a raster by using a polygon shapefile.

The River Architect dataset includes a slackwater zone and zonal statistics help to identify the mean water depth and flow velocity of slackwaters, which are a so-called morphological unit (recall the create-shapefile section).

Background: Instream morphological units aid in describing the geospatial organization of fluvial landforms, which play an important role in ecohydraulic analyses and river restoration. For instance, pool units describe deep water zones with low flow velocity, riffle units are typically characterized by shallow water depths and high velocity, and slackwater units are shallow flow zones with low flow velocity (many juvenile fish love slackwaters). Wyrick and Pasternack (2014) introduce the delineation of morphological units and an open-access summary can be found in the Appendix Sect. 5 in Schwindt et al. (2020).

To analyze a visually apparent slackwater unit, we can draw a polygon in a new shapefile that delineates morphological units. The following figures guide through the creation of a polygon shapefile and the delineation of the riffle with QGIS. Start with opening QGIS and create a new project. Import the water depth and flow velocity rasters showing the slow and shallow water zone. Then follow the workflow in the illustrations below.

Finalize the drawing by clicking on the Save Edits disk button (between Toggle Editing and Add Polygon). Just in case, the slackwater delineation polygon shapefile is also available at the jupyter-python repository of this eBook.

Recall: The new polygon is not saved as long as the edits are not saved. That means: Regularly saving edits when drawing features in QGIS.

Zonal statistics can be calculated using the gdal and ogr libraries, but this is a little cumbersome. The rasterio (conda install -c conda-forge rasterio) library provides a much more convenient method to calculate zonal statistics with its rasterstats.zonal_stats(SHP-FILE, RASTER, STATSTICS-TYPES) method. With zonal_stats, we can easily obtain statistics of the water depth and flow velocity rasters in the limits of the new slackwater polygon.

Recall that both rasters are in the U.S. customary unit system (i.e., feet and feet per second). More statistics can be calculated with zonal_stats:
min, max, mean, count, sum, std, median, majority, minority, unique, range, nodata, percentile_<q> (where <q> can be any float number between 0 and 100).

In addition, user-defined statistics can be added, where the numpy.ma module is particularly useful with its array handling capacities, such as for transposing or specifying statistics along an axis. For instance, we can define a specific function to calculate the standard deviation as follows:

To use the raster_std function in zonal_stats, write something similar to this:

Clip a Raster¶

The above-introduced rasterstats.zonal_stats method works with "Mini-Rasters", which represent clips of the input raster to the user-defined polygon shapefile. Moreover, a mini-raster itself can be obtained by defining the optional keyword argument raster_out=True. In the case that we want to get the original raster clipped without and statistical operation, we can use a little trick by defining an additional statistics function that returns the original array:

With raster_out=True and the original() function, we can retrieve the clipped raster in the following array types:

• mini_raster_array - a clipped and masked numpy array,
• mini_raster_affine - a transformation as an Affine object (i.e., with affine 6-parameter transformation), and
• mini_raster_nodata - NoData values.

The following code block illustrates the usage for retrieving a numpy array:

Tip: Use the above-shown methods to assign a projection and save the clipped array as a GeoTIFF raster. The functions are also implemented in [flusstools] (more precisely, in the script [flusstools.geotools.raster_mgmt.create_raster](https://flusstools.readthedocs.io/en/latest/geotools.html#module-flusstools.geotools.raster_mgmt, which can be used as geotools.raster_mgmt()).

Slope / Aspect Maps and Built-in Command Line Scripts¶

Hillslope maps are an important parameter in hydraulics, hydrology, and ecology. For instance, the slope determines the flow direction of water and it is also a criterion for delineating the habitat of many species. To calculate hill slope gradients and directions, gdal has a command-line tool called gdaldem , which requires a DEM (Digital Elevation Model) raster. The general use of gdaldem in the command line is (arguments in brackets are optional):

gdaldem slope input_dem output_slope_map  [-p use percent slope (default=degrees)] [-s scale* (default=1)] [-alg ZevenbergenThorne] [-compute_edges] [-b Band (default=1)] [-of format] [-co "NAME=VALUE"]* [-q]

To call the command-line tool, we can either use a terminal/command prompt or Python's standard library subprocess. The following code block illustrates the usage of the gdaldem command line tool through subprocess.call() to create a slope raster (in percent) from the River Architect sample data's dem.tif. Note that subprocess.call() returns 0 if the command execution was successful and any other return value indicates an error.

In addition to the absolute slope (i.e., inclination in degrees), or instead of the slope, it can be important to know the slope direction (e.g., inclination toward south, west, east, or north). A raster indicating the slope direction is called aspect raster, where south corresponds to 0° (and 360°), west to 90°, north to 180°, and east to 270°. An aspect raster can also be created with gdaldem:

gdaldem aspect input_dem output_aspect_map [-trigonometric] [-zero_for_flat] [-alg ZevenbergenThorne] [-compute_edges] [-b Band (default=1)] [-of format] [-co "NAME=VALUE"]* [-q]

To create an aspect raster of the River Architect sample data DEM run:

Least Cost Path between Pixels (and another Way of Reprojection)¶

Ecohydraulic Background¶

Least cost paths are important to plan efficient routes for navigation (e.g., in a car) and they can also be helpful in ecohydraulics. Let's take for a moment the position of a fish that, after a flood with decreasing discharge, wants to swim as fast as possible from the floodplain back into the main channel where there is enough water. In the figure below, point 1 shows the starting point on the floodplain and point 2 the destination in the main channel. The reddish background represents the above-produced slope raster (slope-percent.tif) and the water depth at average annual discharge is colored in blue.

Naturally, the path of least costs corresponds to the path of the steepest, monotonously downward-facing slope and we will assume that a fish is able to find it.

A note on hydro-peaking: In many rivers, fish face the daily challenge of escaping from the deadly trap of lateral depressions. The reason for this is that many hydroelectric power plants cause abrupt fluctuations in discharge due to fluctuations in energy demand and production (so-called hydro-peaking). As a result, there is a high stranding risk for fish in many regulated rivers that are subjected to hydro-peaking. The following figure illustrates stranding risk zones as a function of discharge (in cubic feet per second) at the lower Yuba River (California, USA).

Image source: The River Architect Wiki / Kenneth Larrieu

Functions and Libraries Involved¶

The skimage (scikit-image) library (see Other packages in the Open source libraries section) provides with skimage.graph.route_through_array a smart method to calculate the least cost path by summing up pixel-wise connections from point 1 to point 2.

Here is how it works: Assume a numpy array (e.g., with random slope values) that looks like this:

To find the fastest way from array index [0][0] (top left point_1 = (0, 0)) to array index [2][4] (bottom right point point_2 = (2, 4)), we can use route_through_array() to get a list (least_cost_path_indices) with the array coordinates of the path to go, and the costs (weight) involved (sum of all pixels of the least cost path):

To integrate the least cost path list into an array that we can rasterize (geotools.create_raster()), we can paste least_cost_path_indices into a numpy-zeros array of the original slope raster (image) as a transposed list.

In practice, the slope raster is georeferenced, and therefore, we have to use pixel coordinates relative to the coordinate system origin. For this purpose, we need two more functions:

• One function to calculate the pixel-index related offset that we will name coords2offset: The coords2offset() function will return the x-y shift in the form of "number of pixels" (two integers, one for x and one for y shift).
• The above-defined get_srs() function (i.e., geotools.get_srs()).

The coords2offset() function looks like this:

Tip: The coords2offset() function is also available in flusstools with more robust raise-exception notation: use with geotools.coords2offset() and have a look into the script located in flusstools.geotools.dataset_mgmt.

The coords2offset() function converts a raster array (e.g., produced with the above-defined geotools.raster2array() function) into an array that can be used with route_through_array() with the following workflow:

1. Use the raster's geo_transform (gdal.Dataset.GetGeoTransform = (origin_x, pixel_width, 0, origin_y, 0, pixel_height)) and the start and end point coordinates (i.e, start_coord of point 1 and stop_coord of point 2) in coords2offset() to get their pixel indices (start_index_x, start_index_y, stop_index_x, and stop_index_y) in the raster array.
2. Replace np.nan values in the raster array with values that are higher than the maximum value of the array. Do not use zeros, because we want to exclude the np.nan pixels from the least cost path later by overwriting np.nan with very high pixel costs.
3. Use route_through_array() as above explained with the optional arguments geometric=True (use the MCP_Geometric class rather than MCP_base to calculate costs) and fully_connected=True (enables using diagonal pixels as direct neighbors).
4. Integrate the least cost path list (index_path) into a numpy-zeros array (child of raster_array, as explained above) and return the path_array.

Application¶

Recall, we defined the following functions (all are available in flusstools through from flusstools import geotools) that we can use for the calculation of the least cost path to get from point 1 to point 2 in the slope-percent.tif raster:

• geotools.raster2array()
• geotools.create_path_array()
• geotools.get_srs()
• geotools.create_raster()

The below code block uses these functions as follows:

1. Define input (slope-percent.tif) and output (least_cost.tif) raster names (with directories).
2. Define the coordinates of points 1 and 2 as tuples (x, y) in the EPSG:6418 projection.
3. Load the input raster(src_raster), its band as array (raster_array), and geotransformation (geo_transform) with the raster2array() function.
4. Get the least cost path indicated with ones in an array of zeros ( i.e., an on-off path_array) with the create_path_array() function.
5. Get the osgeo.osr.SpatialReference of the input raster (src_raster = osgeo.gdal.Dataset("slope-percent.tif")).
6. Create the least cost path GeoTIFF raster with the create_raster() function as a gdal.GDT_Byte band.

Legitimately, you may wonder whether it was better to represent the least cost path as a line. Of course, that is correct. However, this operation is a conversion of a raster into a line shapefile, which is explained in the next section on geodata conversion. Curious readers can also directly use the raster2line() function of flusstools (or have a look at it in the flusstools.geotools/dataset_mgmt.py script).

Exercise: Familiarize with raster handling in the geospatial ecohydraulics exercise.