pyeo.coordinate_manipulation

Contains a set of functions for transforming spatial coordinates between projections and pixel indicies.

Unless otherwise stated, all functions assume that any geometry, rasters and shapefiles are using the same projection. If they are not, there may be unexpected errors.

Some of these functions call for an AOI shapefile. This is a single-layer shapefile containing only the geometry of one polygon.

These functions mainly work on the objects provided by the ogr and gdal libraries. If you wish to use them in your own processing, a gdal.Image object is usually the output from gdal.Open() and an ogr.Geometry object can be obtained from a well-known text (wkt) string using the snippet object=ogr.ImportFromWkt("mywkt"). For more information on wkt, see https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry and the “QuickWKT” QGIS plugin.

Key functions

pixel_bounds_from_polygon() Gets the indicies of the bounding box of a polygon within a raster

Function reference

pyeo.coordinate_manipulation.align_bounds_to_whole_number(bounding_box)

Creates a new bounding box with it’s height and width rounded to the nearest whole number.

Parameters

bounding_box (ogr.Geometry) – An ogr.Geometry object containing a raster’s bounding box as a polygon.

Returns

bounds_poly – An ogr.Geometry object containing the aligned bounding box.

Return type

ogr.Geometry

pyeo.coordinate_manipulation.check_overlap(raster, aoi)

A test to see if a raster and an AOI overlap. :param raster: A gdal.Image object :type raster: gdal.Image :param aoi: A ogr.Dataset object containing a single polygon :type aoi: ogr.Dataset

Returns

is_overlapped – True if the raster and the polygon overlap, otherwise False.

Return type

bool

pyeo.coordinate_manipulation.floor_to_resolution(input, resolution)

Returns input rounded DOWN to the nearest multiple of resolution. Used to prevent float errors on pixel boarders.

Parameters
  • input (number) – The value to be rounded

  • resolution (number) – The resolution

Returns

output – The largest value between input and 0 that is wholly divisible by resolution.

Return type

number

Notes

Uses the following formula: input-(input%resolution) If resolution is less than 1, then this assumes that the projection is in decmial degrees and will be rounded to 6dp before flooring. However, it is not recommended to use this function in that situation.

pyeo.coordinate_manipulation.get_aoi_bounds(aoi)

Returns a wkbPolygon geometry with the bounding rectangle of a single-layer shapefile

Parameters

aoi (ogr.Dataset) – An ogr.Dataset object containing a single layer.

Returns

bounds_poly – A polygon containing the bounding rectangle

Return type

ogr.Geometry

pyeo.coordinate_manipulation.get_aoi_intersection(raster, aoi)

Returns a wkbPolygon geometry with the intersection of a raster and a shpefile containing an area of interest

Parameters
  • raster (gdal.Image) – A raster containing image data

  • aoi (ogr.DataSource) – A datasource with a single layer and feature

Returns

intersection – An ogr.Geometry object containing a single polygon with the area of intersection

Return type

ogr.Geometry

pyeo.coordinate_manipulation.get_aoi_size(aoi)

Returns the width and height of the bounding box of an aoi.

Parameters

aoi (ogr.Geometry) – A shapefile containing a single layer with a single polygon

Returns

width, height – Width and height of the bounding box

Return type

number

pyeo.coordinate_manipulation.get_combined_polygon(rasters, geometry_mode='intersect')

Returns a polygon containing the combined boundary of each raster in rasters.

Parameters
  • rasters (List of gdal.Image) – A list of raster objects opened with gdal.Open()

  • geometry_mode ({'intersect', 'union'}) – If ‘intersect’, returns the boundary of the area that all rasters cover. If ‘union’, returns the boundary of the area that any raster covers.

Returns

combination – ogr.Geometry() containing a polygon.

Return type

ogr.Geometry

pyeo.coordinate_manipulation.get_local_top_left(raster1, raster2)

Gets the top-left corner of raster1 in the array of raster 2. Assumes both rasters are in the same projection and units.

Parameters
  • raster1 (gdal.Image) – The raster to get the top-left corner of.

  • raster2 (gdal.Image) – The raster that raster1’s top-left corner is over.

Returns

x_pixel, y_pixel – The indicies of the pixel of top-left corner of raster 1 in raster 2.

Return type

int

pyeo.coordinate_manipulation.get_poly_bounding_rect(poly)

Returns a polygon of the bounding rectangle of an input polygon. :param poly: An ogr.Geometry object containing a polygon :type poly: ogr.Geometry

Returns

bounding_rect – An ogr.Geometry object with a four-point polygon representing the bounding rectangle.

Return type

ogr.Geometry

pyeo.coordinate_manipulation.get_poly_intersection(poly1, poly2)

A functional wrapper for ogr.Geometry.Intersection()

Parameters
  • poly1 (ogr.Geometry) – The two geometries to intersect

  • poly2 (ogr.Geometry) – The two geometries to intersect

Returns

intersection – A polygon of the intersection

Return type

ogr.Geometry

pyeo.coordinate_manipulation.get_poly_size(poly)

Returns the width and height of a bounding box of a polygon

Parameters

poly (ogr.Geometry) – A ogr.Geometry object containing the polygon.

Returns

width, height – Width and height of the polygon

Return type

number

pyeo.coordinate_manipulation.get_raster_bounds(raster)

Returns a wkbPolygon geometry with the bounding rectangle of a raster calculated from its geotransform.

Parameters

raster (gdal.Image) – A gdal.Image object

Returns

boundary – An ogr.Geometry object containing a single wkbPolygon with four points defining the bounding rectangle of the raster.

Return type

ogr.Geometry

Notes

Bounding rectangle is obtained from raster.GetGeoTransform(), with the top left corners rounded down to the nearest multiple of of the resolution of the geotransform. This is to avoid rounding errors in reprojected geotransformations.

pyeo.coordinate_manipulation.get_raster_intersection(raster1, raster2)

Returns a wkbPolygon geometry with the intersection of two raster bounding boxes.

Parameters
  • raster1 (gdal.Image) – The gdal.Image

  • raster2 (gdal.Image) – The gdal.Image

Returns

intersections – a ogr.Geometry object containing a single polygon

Return type

ogr.Geometry

pyeo.coordinate_manipulation.get_raster_size(raster)

Return the width and height of a raster, in that raster’s units.

Parameters

raster (gdal.Image) – A gdal.Image object

Returns

width, height – A tuple containing (width, height)

Return type

number

pyeo.coordinate_manipulation.multiple_intersection(polygons)

Takes a list of polygons and returns a geometry representing the intersection of all of them

Parameters

polygons (list of ogr.Geometry) – A list of ogr.Geometry objects, each containing a single polygon.

Returns

intersection – An ogr.Geometry object containing a single polygon

Return type

ogr.Geometry

pyeo.coordinate_manipulation.multiple_union(polygons)

Takes a list of polygons and returns a polygon of the union of their perimeter

Parameters

polygons (list of ogr.Geometry) – A list of ogr.Geometry objects, each containing a single polygon.

Returns

union – An ogr.Geometry object containing a single polygon

Return type

ogr.Geometry

pyeo.coordinate_manipulation.pixel_bounds_from_polygon(raster, polygon)

Returns the bounding box of the overlap between a raster and a polygon in the raster

Parameters
  • raster (gdal.Image) – A gdal raster object

  • polygon (ogr.Geometry or str) – A ogr.Geometry object or wkt string containing a single polygon

Returns

x_min, x_max, y_min, y_max – The bounding box in the pixels of the raster

Return type

int

pyeo.coordinate_manipulation.pixel_to_point_coordinates(pixel, GT)

Given a pixel and a geotransformation, returns the picaltion of that pixel’s top left corner in the projection used by the geotransform. NOTE: At present, this takes input in the form of y,x! This is opposite to the output of point_to_pixel_coordinates!

Parameters
  • pixel (iterable) – A tuple (y, x) of the coordinates of the pixel

  • GT (iterable) – A six-element numpy array containing a geotransform

Returns

Xgeo, Ygeo – The geographic coordinates of the top-left corner of the pixel.

Return type

number

pyeo.coordinate_manipulation.point_to_pixel_coordinates(raster, point, oob_fail=False)
Returns a tuple (x_pixel, y_pixel) in a georaster raster corresponding to the geographic point in a projection.

Assumes raster is north-up non rotated.

Parameters
  • raster (gdal.Image) – A gdal raster object

  • point (str, iterable or ogr.Geometry) –

    One of:

    A well-known text string of a single point An iterable of the form (x,y) An ogr.Geometry object containing a single point

Returns

Return type

A tuple of (x_pixel, y_pixel), containing the indicies of the point in the raster.

Notes

The equation is a rearrangement of the section on affinine geotransform in http://www.gdal.org/gdal_datamodel.html

pyeo.coordinate_manipulation.reproject_geotransform(in_gt, old_proj_wkt, new_proj_wkt)

Reprojects a geotransform from the old projection to a new projection. See [https://gdal.org/user/raster_data_model.html]

Parameters
  • in_gt (array_like) – A six-element numpy array, usually an output from gdal_image.GetGeoTransform()

  • old_proj_wkt (str) – The projection of the old geotransform in well-known text.

  • new_proj_wkt (str) – The projection of the new geotrasform in well-known text.

Returns

out_gt – The geotransform in the new projection

Return type

array_like

pyeo.coordinate_manipulation.reproject_vector(in_path, out_path, dest_srs)

Reprojects a vector file to a new SRS. Simple wrapper for ogr2ogr.

Parameters
  • in_path (str) – Path to the vector file

  • out_path (str) – Path to output vector file

  • dest_srs (osr.SpatialReference or str) – The spatial reference system to reproject to

pyeo.coordinate_manipulation.write_geometry(geometry, out_path, srs_id=4326)

Saves the geometry in an ogr.Geometry object to a shapefile.

Parameters
  • geometry (ogr.Geometry) – An ogr.Geometry object

  • out_path (str) – The location to save the output shapefile

  • srs_id (int or str, optional) – The projection of the output shapefile. Can be an EPSG number or a WKT string.

Notes

The shapefile consists of one layer named ‘geometry’.