GeoTools#

Geospatial Functions for Hydraulics and Morphodynamics

The GeoTools (flusstools.geotools) modules provide Python3 functions for many sorts of river-related analyses with geospatial data (e.g., for working with numerical model input and output). The package is intended as support material for the hydro-informatics eBook.

Usage#

Import#

Import geotools from flusstools:

from flusstools import geotools as geo

Example (code block)#

from flusstools import geotools as geo
raster, array, geo_transform = geo.raster2array("/sample-data/froude.tif")
type(raster)
# >>> <class 'osgeo.gdal.Dataset'>
type(array)
# >>> <class 'numpy.ndarray'>
type(geo_transform)
# >>> <class 'tuple'>
print(geo_transform)
# >>> (6748604.7742, 3.0, 0.0, 2207317.1771, 0.0, -3.0)

Example (showcase)#

A showcase is provided with the ROOT/examples/geotools-showcase/georeference_tifs.py script that illustrates geo-referencing tif images that do not have a projection assigned.

Code structure#

The following diagram highlights function locations in Python scripts and how those are linked to each other.

structure

Fig. 1 Diagram of the code structure (needs to be updated).#

Script and function docs#

Package Head: geotools#

geotools is a package for creating, modifying, and transforming geospatial datasets.

flusstools.geotools.geotools.float2int(raster_file_name, band_number=1)[source]#

Converts a float number raster to an integer raster (required for converting a raster to a polygon shapefile).

Parameters
  • raster_file_name (str) – Target file name, including directory; must end on ".tif".

  • band_number (int) – The raster band number to open (default: 1).

Returns

"path/to/ew_raster_file.tif"

Return type

str

flusstools.geotools.geotools.raster2line(raster_file_name, out_shp_fn, pixel_value, max_distance_method='simplified')[source]#

Converts a raster to a line shapefile, where pixel_value determines line start and end points.

Parameters
  • raster_file_name (str) – of input raster file name, including directory; must end on ".tif".

  • out_shp_fn (str) – of target shapefile name, including directory; must end on ".shp".

  • pixel_value (int or float) – Pixel values to connect.

  • max_distance_method – change to (pixel) "width" or "height" to force lines to exactly follow pixels (no triangulation).

flusstools.geotools.geotools.raster2polygon(file_name, out_shp_fn, band_number=1, field_name='values')[source]#

Converts a raster to a polygon shapefile.

Parameters
  • file_name (str) – Target file name, including directory; must end on ".tif"

  • out_shp_fn (str) – Shapefile name (with directory e.g., "C:/temp/poly.shp")

  • band_number (int) – Raster band number to open (default: 1)

  • field_name (str) – Field name where raster pixel values will be stored (default: "values")

  • add_area – If True, an “area” field will be added, where the area in the shapefiles unit system is calculated (default: False)

flusstools.geotools.geotools.rasterize(in_shp_file_name, out_raster_file_name, pixel_size=10, no_data_value=-9999, rdtype=gdal.GDT_Float32, overwrite=True, interpolate_gap_pixels=False, **kwargs)[source]#

Converts any ESRI shapefile to a raster.

Parameters
  • in_shp_file_name (str) – A shapefile name (with directory e.g., "C:/temp/poly.shp")

  • out_raster_file_name (str) – Target file name, including directory; must end on ".tif"

  • pixel_size (float) – Pixel size as multiple of length units defined in the spatial reference (default: 10)

  • no_data_value (int OR float) – Numeric value for no-data pixels (default: -9999)

  • rdtype (gdal.GDALDataType) – The raster data type (default: gdal.GDT_Float32 (32 bit floating point)

  • overwrite (bool) – Overwrite existing files (default: True)

  • interpolate_gap_pixels (bool) – Fill empty pixels that are not touched by a shapefile element with interpolated values (default: False)

Keyword Arguments
  • field_name (str) – Name of the shapefile’s field with values to burn to raster pixel values.

  • radius1 (float) – Define the x-radius for interpolating pixels (default: -1, corresponding to infinity). Only applicable with interpolate_gap_pixels.

  • radius2 (float) – Define the y-radius for interpolating pixels (default: -1, corresponding to infinity). Only applicable with interpolate_gap_pixels.

  • power (float) – Power of the function for interpolating pixel values (default: 1.0, corresponding to linear).

  • smoothing (float) – Smoothing parameter for interpolating pixel values (default: 0.0).

  • min_points (int) – Minimum number of points to use for interpolation. If the interpolator cannot find at least min_points for a pixel, it assigns a no_data value to that pixel (default: 0).

  • max_points (int) – Maximum number of points to use for interpolation. The interpolator will not use more than max_points closest points to interpolate a pixel value (default: 0).

Hints:

More information on pixel value interpolation: interpolate_gap_pixels=True interpolates values at pixels that are not touched by any las point. The pixel value interpolation uses gdal_grid (i.e., its Python bindings through gdal.Grid()). Control the interpolation parameters with the keyword arguments radius1, radius2, power, max_points, min_points, and smoothing..

Returns

Creates the GeoTIFF raster defined with out_raster_file_name (success: 0, otherwise None).

Return type

int

Raster Management raster_mgmt#

flusstools.geotools.raster_mgmt.clip_raster(polygon, in_raster, out_raster)[source]#

Clips a raster to a polygon.

Parameters
  • polygon (str) – A polygon shapefile name, including directory; must end on ".shp".

  • in_raster (str) – Name of the raster to be clipped, including its directory.

  • out_raster (str) – Name of the target raster, including its directory.

Returns

Creates a new, clipped raster defined with out_raster.

Return type

None

flusstools.geotools.raster_mgmt.create_raster(file_name, raster_array, bands=1, origin=None, epsg=4326, pixel_width=10.0, pixel_height=10.0, nan_val=-9999.0, rdtype=gdal.GDT_Float32, geo_info=False, rotation_angle=None, shear_pixels=True, options=['PROFILE=GeoTIFF'])[source]#

Converts an ndarray (numpy.array) to a GeoTIFF raster.

Parameters
  • file_name (str) – Target file name, including directory; must end on ".tif".

  • raster_array (ndarray or list) – 2d-array (no bands) or list (bands) of 2-darrays of values to rasterize. If a list of 2d-arrays is provided, the length of the list will correspond to the number of bands added to the raster (supersedes bands).

  • bands (int) – Number of bands to write to the raster (default: 1).

  • origin (tuple) – Coordinates (x, y) of the origin.

  • epsg (int) – EPSG:XXXX projection to use (default: 4326).

  • pixel_height (float OR int) – Pixel height as multiple of the base units defined with the EPSG number (default: 10 meters).

  • pixel_width (float OR int) – Pixel width as multiple of the base units defined with the EPSG number (default: 10 meters).

  • nan_val (int or float) – No-data value to be used in the raster. Replaces non-numeric and np.nan in the ndarray. (default: geoconfig.nan_value).

  • rdtypegdal.GDALDataType raster data type (default: gdal.GDT_Float32 (32 bit floating point).

  • geo_info (tuple) – Defines a gdal.DataSet.GetGeoTransform object and supersedes origin, pixel_width, pixel_height (default: False).

  • rotation_angle (float) – Rotate (in degrees) not North-up rasters. The default value (0) corresponds to north-up (only modify if you know what you are doing).

  • shear_pixels (bool) – Use with rotation_angle to shear pixels as well (default: True).

  • options (list) – Raster creation options - default is [‘PROFILE=GeoTIFF’]. Add ‘PHOTOMETRIC=RGB’ to create an RGB image raster.

Returns

0 if successful, otherwise -1.

Return type

int

Hint

For processing airborne imagery, the rotation_angle corresponds to the bearing angle of the aircraft with reference to true, not magnetic North.

flusstools.geotools.raster_mgmt.open_raster(file_name, band_number=1)[source]#

Opens a raster file and accesses its bands.

Parameters
  • file_name (str) – The raster file directory and name.

  • band_number (int) – The Raster band number to open (default: 1).

Returns

A raster dataset a Python object. osgeo.gdal.Band: The defined raster band as Python object.

Return type

osgeo.gdal.Dataset

flusstools.geotools.raster_mgmt.raster2array(file_name, band_number=1)[source]#

Extracts a numpy ndarray from a raster.

Parameters
  • file_name (str) – Target file name, including directory; must end on ".tif".

  • band_number (int) – The raster band number to open (default: 1).

Returns

three-elements of [osgeo.DataSet of the raster, numpy.ndarray of the raster band_numer (input) where no-data values are replaced with np.nan, osgeo.GeoTransform of the original raster]

Return type

list

flusstools.geotools.raster_mgmt.remove_tif(file_name)[source]#

Removes a GeoTIFF and its dependent files (e.g., xml).

Parameters

file_name (str) – Directory (path) and name of a GeoTIFF

Returns

Removes the provided file_name and all dependencies.

Return type

None

flusstools.geotools.raster_mgmt.xy_raster_shift(file_name, x_shift, y_shift, bands=1, rdtype=gdal.GDT_Float32, nan_val=-9999.0, compress=True, options=['PROFILE=GeoTIFF'], compress_config=['COMPRESS=LZW', 'TILED=YES'])[source]#

Creates new geotiff raster with shifts in x and y direction. If enabled compresses it also compresses file.

Args: file_name (string): File path of GeoTiff x_shift (float OR int): Shift origin in x direction. *Check that correct units are used. Example: wgs 84 is in degrees y_shift (float OR int): Shift origin in y direction. *Check that correct units are used. Example: wgs 84 is in degrees bands (int): Number of bands default is 1, however check raster to see how many are required. Example: RGBA=4 rdtype: gdal.GDALDataType raster data type (default: gdal.GDT_Float32 (32 bit floating point). nan_val (int or float): No-data value to be used in the raster. Replaces non-numeric and np.nan in the ndarray. (default: geoconfig.nan_value). compress (Bool): If True creates compressed version of the GeoTiff options (list): Raster creation options - default is [‘PROFILE=GeoTIFF’]. Add ‘PHOTOMETRIC=RGB’ to create an RGB image raster. compress_config: (list) Compress creation options - default is [“COMPRESS=LZW”, “TILED=YES”] LZW=Lempel-Ziv-Welch-Algorithm See gdal.Translate for more options

Returns

0 if successful, otherwise -1.

Return type

int

Hint: For drone rasters try bands=4 (rgba) rdtype=gdal.GDT_Byte nan_val=0 options=[‘PROFILE=GeoTIFF’,’PHOTOMETRIC=RGB’])

Bugs: Issues displaying logging

Shapefile Management shp_mgmt#

flusstools.geotools.shp_mgmt.create_shp(shp_file_dir, overwrite=True, *args, **kwargs)[source]#

Creates a new shapefile with an optionally defined geometry type.

Parameters
  • shp_file_dir (str) – of the (relative) shapefile directory (ends on ".shp").

  • overwrite (bool) – If True (default), existing files are overwritten.

  • layer_name (str) – The layer name to be created. If None: no layer will be created.

  • layer_type (str) – Either "point", "line", or "polygon" of the layer_name. If None: no layer will be created.

Returns

An ogr shapefile (Python object)

Return type

osgeo.ogr.DataSource

Hint

Use the layer_name and layer_type kwargs along with each other. Keeping these parameters default is deprecated.

flusstools.geotools.shp_mgmt.get_geom_description(layer)[source]#

Gets the WKB Geometry Type as string from a shapefile layer.

Parameters

layer (osgeo.ogr.Layer) – A shapefile layer.

Returns

WKB (binary) geometry type

Return type

str

flusstools.geotools.shp_mgmt.get_geom_simplified(layer)[source]#
Gets a simplified geometry description (either point, line, or polygon)

as a function of the WKB Geometry Type of a shapefile layer.

Parameters

layer (osgeo.ogr.Layer) – A shapefile layer.

Returns

Either WKT-formatted point, line, or polygon (or unknown if invalid layer).

Return type

str

flusstools.geotools.shp_mgmt.polygon_from_shapepoints(shapepoints, polygon, alpha=nan)[source]#

Creates a polygon around a cloud of shapepoints.

Parameters
  • shapepoints (str) – Point shapefile name, including its directory.

  • polygon (str) – Target shapefile filename, including its directory.

  • alpha (float) – Coefficient to adjust; the lower it is, the more slim will be the polygon.

Returns

Creates the polygon shapefile defined with polygon.

Return type

None

flusstools.geotools.shp_mgmt.verify_shp_name(shp_file_name, shorten_to=13)[source]#

Ensure that the shapefile name does not exceed 13 characters. Otherwise, the function shortens the shp_file_name length to shorten_to=N characters.

Parameters
  • shp_file_name (str) – A shapefile name (with directory e.g., "C:/temp/poly.shp").

  • shorten_to (int) – The number of characters the shapefile name should have (default: 13).

Returns

A shapefile name (including path if provided) with a length of shorten_to.

Return type

str

Projection Management srs_mgmt#

flusstools.geotools.srs_mgmt.get_esriwkt(epsg)[source]#

Gets esriwkt-formatted spatial references with epsg code online.

Parameters

epsg (int) – EPSG Authority Code

Returns

An esriwkt string (if an error occur, the default epsg=``4326`` is used).

Return type

str

Example

Call this function with get_esriwkt(4326) to get a return, such as 'GEOGCS["GCS_WGS_1984",DATUM[...],...].

Hint

This function requires an internet connection: Loads spatial reference codes as "https://spatialreference.org/ref/sr-org/{0}/esriwkt/".format(epsg) For instance, epsg=3857 yields "https://spatialreference.org/ref/sr-org/3857/esriwkt/"

flusstools.geotools.srs_mgmt.get_srs(dataset)[source]#

Gets the spatial reference of any gdal.Dataset.

Parameters

dataset (gdal.Dataset) – A shapefile or raster.

Returns

A spatial reference object.

Return type

osr.SpatialReference

flusstools.geotools.srs_mgmt.get_wkt(epsg, wkt_format='esriwkt')[source]#

Gets WKT-formatted projection information for an epsg code using the osr library.

Parameters
  • epsg (int) – epsg Authority code

  • wkt_format (str) – of wkt format (default is esriwkt for shapefile projections)

Returns

WKT (if error: returns default corresponding to epsg=4326).

Return type

str

flusstools.geotools.srs_mgmt.make_prj(shp_file_name, epsg)[source]#

Generates a projection file for a shapefile.

Parameters
  • shp_file_name (str) – of a shapefile name (with directory e.g., "C:/temp/poly.shp").

  • epsg (int) – EPSG Authority Code

Returns

Creates a projection file (.prj) in the same directory and with the same name of shp_file_name.

Return type

None

flusstools.geotools.srs_mgmt.reproject(source_dataset, new_projection_dataset)[source]#

Re-projects a dataset (raster or shapefile) onto the spatial reference system of a (shapefile or raster) layer.

Parameters
  • source_dataset (gdal.Dataset) – Shapefile or raster.

  • new_projection_dataset (gdal.Dataset) – Shapefile or raster with new projection info.

Returns

If the source is a raster, the function creates a GeoTIFF in same directory as source_dataset with a "_reprojected" suffix in the file name. If the source is a shapefile, the function creates a shapefile in same directory as source_dataset with a "_reprojected" suffix in the file name.

Return type

None

flusstools.geotools.srs_mgmt.reproject_raster(source_dataset, source_srs, target_srs)[source]#

Re-projects a raster dataset. This function is called by the reproject function.

Parameters
  • source_dataset (osgeo.ogr.DataSource) – Instantiates with an ogr.Open(SHP-FILE).

  • source_srs (osgeo.osr.SpatialReference) – Instantiates with get_srs(source_dataset)

  • target_srs (osgeo.osr.SpatialReference) – Instantiates with get_srs(DATASET-WITH-TARGET-PROJECTION).

Returns

Creates a new GeoTIFF raster in the same directory where source_dataset lives.

Return type

None

flusstools.geotools.srs_mgmt.reproject_shapefile(source_dataset, source_layer, source_srs, target_srs)[source]#

Re-projects a shapefile dataset. This function is called by the reproject function.

Parameters
  • source_dataset (osgeo.ogr.DataSource) – Instantiates with ogr.Open(SHP-FILE).

  • source_layer (osgeo.ogr.Layer) – Instantiates with source_dataset.GetLayer().

  • source_srs (osgeo.osr.SpatialReference) – Instantiates with get_srs(source_dataset).

  • target_srs (osgeo.osr.SpatialReference) – Instantiates with get_srs(DATASET-WITH-TARGET-PROJECTION).

Returns

Creates a new shapefile in the same directory where source_dataset lives.

Return type

None

Dataset Management (dataset_mgmt)#

flusstools.geotools.dataset_mgmt.coords2offset(geo_transform, x_coord, y_coord)[source]#

Returns x-y pixel offset (inverse of the offset2coords function).

Parameters
  • geo_transform – osgeo.gdal.Dataset.GetGeoTransform() object

  • x_coord (float) – x-coordinate

  • y_coord (float) – y-coordinate

Returns

Number of pixels (offset_x, offset_y), both int .

Return type

tuple

flusstools.geotools.dataset_mgmt.get_layer(dataset, band_number=1)[source]#

Gets a layer=band (RasterDataSet) or layer=ogr.Dataset.Layer of any dataset.

Parameters
  • dataset (osgeo.gdal.Dataset or osgeo.ogr.DataSource) – Either a raster or a shapefile.

  • band_number (int) – Only use with rasters to define a band number to open (default is 1 ).

Returns

{"type": raster or vector or "None", layer": if raster: raster_band, if vector: GetLayer(), else: None}

Return type

dict

flusstools.geotools.dataset_mgmt.offset2coords(geo_transform, offset_x, offset_y)[source]#

Returns x-y coordinates from pixel offset (inverse of coords2offset function).

Parameters
  • geo_transform (osgeo.gdal.Dataset.GetGeoTransform) – The geo transformation to use.

  • offset_x (int) – x number of pixels.

  • offset_y (int) – y number of pixels.

Returns

Two float numbers of x-y-coordinates (x_coord, y_coord).

Return type

tuple

flusstools.geotools.dataset_mgmt.verify_dataset(dataset)[source]#

Verifies if a dataset contains raster or vector data.

Parameters

dataset (osgeo.gdal.Dataset or osgeo.ogr.DataSource) – Dataset to verify.

Returns

Either “unknown”, “raster”, or “vector”.

Return type

str

KML/KML File Management#

Modified script (original: Linwood Creekmore III)

Examples

output to geopandas dataframe (gdf): gdf = kmx2other("my-places.kmz", output="gpd")

plot the new gdf (use %matplotlib inline in notebooks) gdf.plot()

convert a kml-file to a shapefile success = kmx2other("my-places.kml", output="shp")

flusstools.geotools.kml.kmx2other(file, output='df')[source]#

Converts a Keyhole Markup Language Zipped (KMZ) or KML file to a pandas dataframe, geopandas geodataframe, csv, geojson, or ESRI shapefile.

Parameters
  • file (str) – The path to a KMZ or KML file.

  • output (str) – Defines the output type. Valid options are: "shapefile", "shp", "shapefile", or "ESRI Shapefile".

Returns

Success message (use print(kmx2other(...)) to see what the function did.)

Return type

str

Original Classes written by Linwood Creekmore III (modified for flusstools)

Flavored with code blocks from:

class flusstools.geotools.kmx_parser.ModHTMLParser[source]#

A child of HTMLParser, tailored (modified) for kml/kmy parsing.

handle_data(data)[source]#

Generates mapping and series if in_table is True.

Parameters

data (str) – Text lines of data divided by colons.

Returns

Assigns ModHTMLParser.mapping and ModHTMLParser.series attributes

Return type

None

handle_starttag(tag, attrs)[source]#

Enables a table if a table-tag is provided.

Parameters
  • tag (str) – Set to “table” for enabling usage of a table.

  • attrs (list) – List of additional attributes (currently unused).

Returns

Verifies if the tag argument contains the string "table"

Return type

None

class flusstools.geotools.kmx_parser.PlacemarkHandler[source]#

Child of xml.sax.handler.ContentHandler, tailored for handling kml files.

characters(data)[source]#

Adds a line of data to the read-buffer.

Parameters

data (str) –

Returns

None

end_element(name)[source]#

Sets the end (last) element.

Parameters

name (str) –

Returns

None

htmlizer()[source]#

Creates an html file.

Parameters

row (pandas.df) – List of strings for conversion

Returns

An instance of the ModHTMLParser() class

Return type

htmlparser.series

spatializer()[source]#

Converts string objects to spatial Python objects.

Parameters

row (pandas.df) – List of strings for conversion

Returns

None

start_element(name)[source]#

Looks for the first Placemark element in a kml file.

Parameters

name (str) – Name-tag of the element

Returns

None

Shortest Path Finder#

This module is inspired by Michael Diener - read more at

mdiener21/python-geospatial-analysis-cookbook

Example use: create_shortest_path(shp_file_name, start_node_id, end_node_id)

flusstools.geotools.shortest_path.create_shortest_path(line_shp_name, start_node_id, end_node_id)[source]#

Calculates the shortest path from a network of lines.

Parameters
  • line_shp_name (str) – Input shapefile name

  • start_node_id (int) – Start node ID

  • end_node_id (int) – End node ID

Returns

Creates a graph of nodes (coordinate pairs) connecting a start node with an end node in the defined line_shp_name.

Return type

None

flusstools.geotools.shortest_path.get_full_path(path, nx_list_subgraph)[source]#

Creates a numpy array of the line result.

Parameters
  • path (str) – Result of nx.shortest_path

  • nx_list_subgraph (list) – See create_shortest path function

Returns

Coordinate pairs along a path.

Return type

ndarray

flusstools.geotools.shortest_path.get_path(n0, n1, nx_list_subgraph)[source]#

Get path between nodes n0 and n1.

Parameters
  • n0 (int) – Node 1

  • n1 (int) – Node 2

  • nx_list_subgraph (list) – (see create shortest path)

Returns

An array of point coordinates along the line linking these two nodes.

Return type

ndarray

flusstools.geotools.shortest_path.write_geojson(outfilename, indata)[source]#

Creates a new GeoJSON file

Args>

outfilename (str): Name for the output file indata (array): Array to write tyo the geojson file.

Returns

Creates a new GeoJSON file.