GeoTools

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.

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).
rdtype – gdal.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".

Hint

The core function is taken from http://programmingadvent.blogspot.com/2013/06/kmzkml-file-parsing-with-python.html

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.

GeoTools

Contents

GeoTools#

Usage#

Import#

Example (code block)#

Example (showcase)#

Code structure#

Script and function docs#

Package Head: geotools#

Raster Management raster_mgmt#

Shapefile Management shp_mgmt#

Projection Management srs_mgmt#

Dataset Management (dataset_mgmt)#

KML/KML File Management#

Shortest Path Finder#

Package Head: `geotools`#

Raster Management `raster_mgmt`#

Shapefile Management `shp_mgmt`#

Projection Management `srs_mgmt`#

Dataset Management (`dataset_mgmt`)#