auscophub.geomutils¶
Some utility functions for handling the geometries given with Sentinel files.
ESA have done a fairly minimal job with their footprint geometries. The routines in this file are mostly aimed at providing a bit more support beyond that, mostly for things like crossing the International Date Line.
It is assumed that we have access to the GDAL/OGR library, and numpy.
Note on efficiency.¶
These routines all operate on the basic ogr.Geometry object. However, since this is a relatively opaque structure, mostly they will convert to JSON strings and use those to get hold of the coordinates. This means there is a fair bit of converting back and forth, which means we are not at maximum efficiency. There is also quite a bit of changing of projection. However, in the context of AusCopHub, this is all somewhat negligible compared to doing things like reading and copying ESA zip files.
-
auscophub.geomutils.copyGeom(geom)¶ Return a copy of the geometry. OGR does not provide a good method for doing this.
-
auscophub.geomutils.crossesDateline(geom, preferredEpsg)¶ Given a Polygon Geometry object, in lat/long, detect whether it crosses the dateline. Do this in the projection of the preferred EPSG, so we remove (reduce?) the ambiguity about inside/outside.
-
auscophub.geomutils.findCentroid(geom, preferredEpsg)¶ Given a geometry as a lat/long polygon, find the lat/long centroid, by first projecting into the preferred EPSG, so as to avoid discontinuities. The preferredEpsg is one in which the polygon ought to make sense (as found, hopefully, by the findSensibleProjection() function).
Returns a pair [centroidX, centroidY] in lat/long
-
auscophub.geomutils.findSensibleProjection(geom)¶ Given a polygon Geometry object in lat/long, work out what would be a suitable projection to use with this area, in order to avoid things like the international date line wrap-around, or the north/sourth pole discontinuities.
This only makes sense for tiled products, as opposed to long strips which cross multiple zones, etc.
Main possible options are UTM in a suitable zone, UPS when near the poles.
Return the EPSG number of the projection.
-
auscophub.geomutils.geomFromInteriorPoints(coords)¶ The given list of pairs (or 2-d numpy array) is the (x, y) coords of a set of internal points inside a polygon. Returns a MultiPoint Geometry.
-
auscophub.geomutils.geomFromOutlineCoords(coords)¶ The given list of pairs (or 2-d numpy array) is the (x, y) coords of the polygon outline. Return a Polygon ogr.Geometry object.
-
auscophub.geomutils.getCoords(geom)¶ Return the coordinates of the given OGR geometry. Assumes that this is a single polygon, and returns a numpy array of the x, y coords, of shape (numPts, 2).
If the polygon has holes, they will be discarded - this is just the outer polygon.
If the geometry is a MultiPoint geom, also return a 2-d array of coords.
-
auscophub.geomutils.makeTransformations(epsg1, epsg2)¶ Make a pair of ogr.CoordinateTransformation objects, for transforming between the two given EPSG projections.
- Return a tuple
- (tr1to2, tr2to1)
-
auscophub.geomutils.polygonFromInteriorPoints(geom, preferredEpsg)¶ Given a MultiPoint geometry object in lat/long, create a polygon of the convex hull of these points.
First project the lat/long points into the preferred EPSG, so that when we find the convex hull, we are not crossing any discontinuities such as the international date line.
Return a single polygon geometry in lat/long.
-
auscophub.geomutils.preventGdal3axisSwap(sr)¶ Guard the given spatial reference object against axis swapping, when running with GDAL 3. Does nothing if GDAL < 3. Modifies the object in place. See https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order for details of this issue.
-
auscophub.geomutils.splitAtDateline(geom, preferredEpsg)¶ Given a Polygon Geometry object in lat/long, determine whether it crosses the date line, and if so, split it into a multipolygon with a part on either side.
Use the given preferred EPSG to perform calculations.
Return a new Geometry in lat/long.