Repository: JuliaGeo/LibGEOS.jl
Pull request #143 (2022-11-06 20:40:15)

change default context of operations from global context to operand context
This PR makes it more safe and convenient to work with contexts other than the global one. 

Before this PR all operations are by default done in the global context.
For objects that were created using another context, this is (almost?) never what you want.

After this PR all operations are by default done in the context of the operands.
I think this almost always what you want. If the contexts of operands mismatch an error is thrown. 

* This PR does not change any behavior if you always work with the default global context.
* This PR also does not change any behavior if you always pass contexts explicitly.


Patch:
diff --git a/src/geos_operations.jl b/src/geos_operations.jl
index 1ff7e93..8459ba7 100644
--- a/src/geos_operations.jl
+++ b/src/geos_operations.jl
@@ -9,11 +9,11 @@ const GEOMTYPE = Dict{GEOSGeomTypes,Symbol}(
     GEOS_GEOMETRYCOLLECTION => :GeometryCollection,
 )
 
-writegeom(obj::Geometry, wktwriter::WKTWriter, context::GEOSContext = get_global_context()) =
+writegeom(obj::Geometry, wktwriter::WKTWriter, context::GEOSContext = get_context(obj)) =
     _writegeom(obj.ptr, wktwriter, context)
-writegeom(obj::Geometry, wkbwriter::WKBWriter, context::GEOSContext = get_global_context()) =
+writegeom(obj::Geometry, wkbwriter::WKBWriter, context::GEOSContext = get_context(obj)) =
     _writegeom(obj.ptr, wkbwriter, context)
-writegeom(obj::Geometry, context::GEOSContext = get_global_context()) = 
+writegeom(obj::Geometry, context::GEOSContext = get_context(obj)) = 
     _writegeom(obj.ptr, context)
 
 function geomFromGEOS(ptr::GEOSGeom, context::GEOSContext = get_global_context())
@@ -51,20 +51,20 @@ readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = get_global_context())
 # -----
 # Linear referencing functions -- there are more, but these are probably sufficient for most purposes
 # -----
-project(line::LineString, point::Point, context::GEOSContext = get_global_context()) =
+project(line::LineString, point::Point, context::GEOSContext = get_context(line)) =
     project(line.ptr, point.ptr, context)
-projectNormalized(line::LineString, point::Point, context::GEOSContext = get_global_context()) =
+projectNormalized(line::LineString, point::Point, context::GEOSContext = get_context(line)) =
     projectNormalized(line.ptr, point.ptr, context)
-interpolate(line::LineString, dist::Real, context::GEOSContext = get_global_context()) =
+interpolate(line::LineString, dist::Real, context::GEOSContext = get_context(line)) =
     Point(interpolate(line.ptr, dist, context), context)
-interpolateNormalized(line::LineString, dist::Real, context::GEOSContext = get_global_context()) =
+interpolateNormalized(line::LineString, dist::Real, context::GEOSContext = get_context(line)) =
     Point(interpolateNormalized(line.ptr, dist, context), context)
 
 # # -----
 # # Topology operations
 # # -----
 
-buffer(obj::Geometry, dist::Real, quadsegs::Integer = 8, context::GEOSContext = get_global_context()) =
+buffer(obj::Geometry, dist::Real, quadsegs::Integer = 8, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(buffer(obj.ptr, dist, quadsegs, context), context)
 bufferWithStyle(
     obj::Geometry,
@@ -73,35 +73,35 @@ bufferWithStyle(
     endCapStyle::GEOSBufCapStyles = GEOSBUF_CAP_ROUND,
     joinStyle::GEOSBufJoinStyles = GEOSBUF_JOIN_ROUND,
     mitreLimit::Real = 5.0,
-    context::GEOSContext = get_global_context(),) =
+    context::GEOSContext = get_context(obj),) =
     geomFromGEOS(bufferWithStyle(obj.ptr, dist, quadsegs, endCapStyle, joinStyle, mitreLimit, context), context)
 
-envelope(obj::Geometry, context::GEOSContext = get_global_context()) =
+envelope(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(envelope(obj.ptr, context), context)
-envelope(obj::PreparedGeometry, context::GEOSContext = get_global_context()) =
+envelope(obj::PreparedGeometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(envelope(obj.ownedby.ptr, context), context)
-minimumRotatedRectangle(obj::Geometry, context::GEOSContext = get_global_context()) =
+minimumRotatedRectangle(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(minimumRotatedRectangle(obj.ptr, context), context)
-convexhull(obj::Geometry, context::GEOSContext = get_global_context()) =
+convexhull(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(convexhull(obj.ptr, context), context)
-boundary(obj::Geometry, context::GEOSContext = get_global_context()) =
+boundary(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(boundary(obj.ptr, context), context)
-unaryUnion(obj::Geometry, context::GEOSContext = get_global_context()) =
+unaryUnion(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(unaryUnion(obj.ptr, context), context)
-pointOnSurface(obj::Geometry, context::GEOSContext = get_global_context()) =
+pointOnSurface(obj::Geometry, context::GEOSContext = get_context(obj)) =
     Point(pointOnSurface(obj.ptr, context), context)
-centroid(obj::Geometry, context::GEOSContext = get_global_context()) =
+centroid(obj::Geometry, context::GEOSContext = get_context(obj)) =
     Point(centroid(obj.ptr, context), context)
-node(obj::Geometry, context::GEOSContext = get_global_context()) =
+node(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(node(obj.ptr, context), context)
 
-intersection(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+intersection(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     geomFromGEOS(intersection(obj1.ptr, obj2.ptr, context), context)
-difference(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+difference(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     geomFromGEOS(difference(obj1.ptr, obj2.ptr, context), context)
-symmetricDifference(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+symmetricDifference(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     geomFromGEOS(symmetricDifference(obj1.ptr, obj2.ptr, context), context)
-union(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+union(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     geomFromGEOS(union(obj1.ptr, obj2.ptr, context), context)
 
 # # all arguments remain ownership of the caller (both Geometries and pointers)
@@ -123,80 +123,84 @@ union(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context(
 #     result
 # end
 
-simplify(obj::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
+simplify(obj::Geometry, tol::Real, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(simplify(obj.ptr, tol, context), context)
-topologyPreserveSimplify(obj::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
+topologyPreserveSimplify(obj::Geometry, tol::Real, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(topologyPreserveSimplify(obj.ptr, tol, context), context)
-uniquePoints(obj::Geometry, context::GEOSContext = get_global_context()) =
+uniquePoints(obj::Geometry, context::GEOSContext = get_context(obj)) =
     MultiPoint(uniquePoints(obj.ptr, context), context)
-delaunayTriangulationEdges(obj::Geometry, tol::Real = 0.0, context::GEOSContext = get_global_context()) =
+delaunayTriangulationEdges(obj::Geometry, tol::Real = 0.0, context::GEOSContext = get_context(obj)) =
     MultiLineString(delaunayTriangulation(obj.ptr, tol, true, context), context)
-delaunayTriangulation(obj::Geometry, tol::Real = 0.0, context::GEOSContext = get_global_context()) =
+delaunayTriangulation(obj::Geometry, tol::Real = 0.0, context::GEOSContext = get_context(obj)) =
     GeometryCollection(delaunayTriangulation(obj.ptr, tol, false, context), context)
-constrainedDelaunayTriangulation(obj::Geometry, context::GEOSContext = get_global_context()) =
+constrainedDelaunayTriangulation(obj::Geometry, context::GEOSContext = get_context(obj)) =
     GeometryCollection(constrainedDelaunayTriangulation(obj.ptr, context), context)
 
 
-sharedPaths(obj1::LineString, obj2::LineString, context::GEOSContext = get_global_context()) =
+sharedPaths(obj1::LineString, obj2::LineString, context::GEOSContext = get_context(obj1,obj2)) =
     GeometryCollection(sharedPaths(obj1.ptr, obj2.ptr, context), context)
 
 # # Snap first geometry on to second with given tolerance
-snap(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
+snap(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = get_context(obj1,obj2)) =
     geomFromGEOS(snap(obj1.ptr, obj2.ptr, tol, context), context)
 
 # -----
 # Binary predicates
 # -----
 
-disjoint(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+disjoint(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     disjoint(obj1.ptr, obj2.ptr, context)
-touches(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+touches(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     touches(obj1.ptr, obj2.ptr, context)
-intersects(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+intersects(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     intersects(obj1.ptr, obj2.ptr, context)
-crosses(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+crosses(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     crosses(obj1.ptr, obj2.ptr, context)
-within(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+within(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     within(obj1.ptr, obj2.ptr, context)
-Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     Base.contains(obj1.ptr, obj2.ptr, context)
-overlaps(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+overlaps(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     overlaps(obj1.ptr, obj2.ptr, context)
-equals(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
-    equals(obj1.ptr, obj2.ptr, context)
-equalsexact(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
+equalsexact(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = get_context(obj1,obj2)) =
     equalsexact(obj1.ptr, obj2.ptr, tol, context)
-covers(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+covers(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     covers(obj1.ptr, obj2.ptr, context)
-coveredby(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+coveredby(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     coveredby(obj1.ptr, obj2.ptr, context)
 
+function equals(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1))
+    # We don't check context here, as it is probably safe to compare objects from 
+    # different contexts
+    equals(obj1.ptr, obj2.ptr, context)
+end
+
 
 # # -----
 # # Prepared Geometry Binary predicates
 # # -----
 
-prepareGeom(obj::Geometry, context::GEOSContext = get_global_context()) =
+prepareGeom(obj::Geometry, context::GEOSContext = get_context(obj)) =
     PreparedGeometry(prepareGeom(obj.ptr, context), obj)
-Base.contains(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+Base.contains(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepcontains(obj1.ptr, obj2.ptr, context)
-containsproperly(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+containsproperly(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepcontainsproperly(obj1.ptr, obj2.ptr, context)
-coveredby(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+coveredby(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepcoveredby(obj1.ptr, obj2.ptr, context)
-covers(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context())= 
+covers(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2))= 
     prepcovers(obj1.ptr, obj2.ptr, context)
-crosses(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+crosses(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepcrosses(obj1.ptr, obj2.ptr, context)
-disjoint(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+disjoint(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepdisjoint(obj1.ptr, obj2.ptr, context)
-intersects(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+intersects(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepintersects(obj1.ptr, obj2.ptr, context)
-overlaps(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+overlaps(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepoverlaps(obj1.ptr, obj2.ptr, context)
-touches(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+touches(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     preptouches(obj1.ptr, obj2.ptr, context)
-within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     prepwithin(obj1.ptr, obj2.ptr, context)
 
 # # -----
@@ -212,20 +216,20 @@ within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global
 # # -----
 # # Unary predicate - return 2 on exception, 1 on true, 0 on false
 # # -----
-isEmpty(obj::Geometry, context::GEOSContext = get_global_context()) =
+isEmpty(obj::Geometry, context::GEOSContext = get_context(obj)) =
     isEmpty(obj.ptr, context)
-isEmpty(obj::PreparedGeometry, context::GEOSContext = get_global_context()) =
+isEmpty(obj::PreparedGeometry, context::GEOSContext = get_context(obj)) =
     isEmpty(obj.ownedby.ptr, context)
-isSimple(obj::Geometry, context::GEOSContext = get_global_context()) =
+isSimple(obj::Geometry, context::GEOSContext = get_context(obj)) =
     isSimple(obj.ptr, context)
-isRing(obj::Geometry, context::GEOSContext = get_global_context()) =
+isRing(obj::Geometry, context::GEOSContext = get_context(obj)) =
     isRing(obj.ptr, context)
-isValid(obj::Geometry, context::GEOSContext = get_global_context()) =
+isValid(obj::Geometry, context::GEOSContext = get_context(obj)) =
     isValid(obj.ptr, context)
-hasZ(obj::Geometry, context::GEOSContext = get_global_context()) =
+hasZ(obj::Geometry, context::GEOSContext = get_context(obj)) =
     hasZ(obj.ptr, context)
 
-isClosed(obj::LineString, context::GEOSContext = get_global_context()) =
+isClosed(obj::LineString, context::GEOSContext = get_context(obj)) =
     isClosed(obj.ptr, context) # Call only on LINESTRING
 
 # # -----
@@ -259,7 +263,7 @@ isClosed(obj::LineString, context::GEOSContext = get_global_context()) =
 # # -----
 
 # Gets the number of sub-geometries
-numGeometries(obj::Geometry, context::GEOSContext = get_global_context()) =
+numGeometries(obj::Geometry, context::GEOSContext = get_context(obj)) =
     numGeometries(obj.ptr, context)
 
 # # Call only on GEOMETRYCOLLECTION or MULTI*
@@ -277,23 +281,23 @@ numGeometries(obj::Geometry, context::GEOSContext = get_global_context()) =
 # end
 # getGeometries(ptr::GEOSGeom) = GEOSGeom[getGeometry(ptr, i) for i=1:numGeometries(ptr)]
 # Gets sub-geomtry at index n or a vector of all sub-geometries
-getGeometry(obj::Geometry, n::Integer, context::GEOSContext = get_global_context()) =
+getGeometry(obj::Geometry, n::Integer, context::GEOSContext = get_context(obj)) =
     geomFromGEOS(getGeometry(obj.ptr, n, context), context)
-getGeometries(obj::Geometry, context::GEOSContext = get_global_context()) =
+getGeometries(obj::Geometry, context::GEOSContext = get_context(obj)) =
     [geomFromGEOS(gptr, context) for gptr in getGeometries(obj.ptr, context)]
 
 # Converts Geometry to normal form (or canonical form).
-normalize!(obj::Geometry, context::GEOSContext = get_global_context()) =
+normalize!(obj::Geometry, context::GEOSContext = get_context(obj)) =
     normalize!(obj.ptr, context)
 
 # LinearRings in Polygons
-numInteriorRings(obj::Polygon, context::GEOSContext = get_global_context()) =
+numInteriorRings(obj::Polygon, context::GEOSContext = get_context(obj)) =
     numInteriorRings(obj.ptr, context)
-interiorRing(obj::Polygon, n::Integer, context::GEOSContext = get_global_context()) =
+interiorRing(obj::Polygon, n::Integer, context::GEOSContext = get_context(obj)) =
     LinearRing(interiorRing(obj.ptr, n, context), context)
-interiorRings(obj::Polygon, context::GEOSContext = get_global_context()) =
+interiorRings(obj::Polygon, context::GEOSContext = get_context(obj)) =
     map(LinearRing, interiorRings(obj.ptr, context))
-exteriorRing(obj::Polygon, context::GEOSContext = get_global_context()) =
+exteriorRing(obj::Polygon, context::GEOSContext = get_context(obj)) =
     LinearRing(exteriorRing(obj.ptr, context), context)
 
 # # Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
@@ -321,33 +325,33 @@ exteriorRing(obj::Polygon, context::GEOSContext = get_global_context()) =
 #     result
 # end
 
-numPoints(obj::LineString, context::GEOSContext = get_global_context()) =
+numPoints(obj::LineString, context::GEOSContext = get_context(obj)) =
     numPoints(obj.ptr, context) # Call only on LINESTRING
-startPoint(obj::LineString, context::GEOSContext = get_global_context()) =
+startPoint(obj::LineString, context::GEOSContext = get_context(obj)) =
     Point(startPoint(obj.ptr, context), context) # Call only on LINESTRING
-endPoint(obj::LineString, context::GEOSContext = get_global_context()) =
+endPoint(obj::LineString, context::GEOSContext = get_context(obj)) =
     Point(endPoint(obj.ptr, context), context) # Call only on LINESTRING
 
 # # -----
 # # Misc functions
 # # -----
 
-area(obj::Geometry, context::GEOSContext = get_global_context()) =
+area(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomArea(obj.ptr, context)
-geomLength(obj::Geometry, context::GEOSContext = get_global_context()) =
+geomLength(obj::Geometry, context::GEOSContext = get_context(obj)) =
     geomLength(obj.ptr, context)
 
-distance(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
+distance(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) =
     geomDistance(obj1.ptr, obj2.ptr, context)
-hausdorffdistance(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) = 
+hausdorffdistance(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2)) = 
     hausdorffdistance(obj1.ptr, obj2.ptr, context)
 
-hausdorffdistance(obj1::Geometry, obj2::Geometry, densify::Real,context::GEOSContext = get_global_context()) =
+hausdorffdistance(obj1::Geometry, obj2::Geometry, densify::Real,context::GEOSContext = get_context(obj1,obj2)) =
     hausdorffdistance(obj1.ptr, obj2.ptr, densify, context)
 
 # Returns the closest points of the two geometries.
 # The first point comes from g1 geometry and the second point comes from g2.
-function nearestPoints(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context())
+function nearestPoints(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_context(obj1,obj2))
     points = nearestPoints(obj1.ptr, obj2.ptr, context)
     if points == C_NULL
         return Point[]
@@ -361,22 +365,22 @@ end
 # # Precision functions
 # # -----
 
-getPrecision(obj::Geometry, context::GEOSContext = get_global_context()) =
+getPrecision(obj::Geometry, context::GEOSContext = get_context(obj)) =
     getPrecision(obj.ptr, context)
 setPrecision(
     obj::Geometry,
     grid::Real;
     flags = GEOS_PREC_VALID_OUTPUT,
-    context::GEOSContext = get_global_context(),) =
+    context::GEOSContext = get_context(obj),) =
     setPrecision(obj.ptr, grid::Real, flags, context)
 
 # ----
 #  Geometry information functions
 # ----
 
-getXMin(obj::Geometry, context::GEOSContext = get_global_context()) = getXMin(obj.ptr, context)
-getYMin(obj::Geometry, context::GEOSContext = get_global_context()) = getYMin(obj.ptr, context)
-getXMax(obj::Geometry, context::GEOSContext = get_global_context()) = getXMax(obj.ptr, context)
-getYMax(obj::Geometry, context::GEOSContext = get_global_context()) = getYMax(obj.ptr, context)
+getXMin(obj::Geometry, context::GEOSContext = get_context(obj)) = getXMin(obj.ptr, context)
+getYMin(obj::Geometry, context::GEOSContext = get_context(obj)) = getYMin(obj.ptr, context)
+getXMax(obj::Geometry, context::GEOSContext = get_context(obj)) = getXMax(obj.ptr, context)
+getYMax(obj::Geometry, context::GEOSContext = get_context(obj)) = getYMax(obj.ptr, context)
 
 # TODO 02/2022: wait for libgeos release beyond 3.10.2 which will in include GEOSGeom_getExtent_r
diff --git a/src/geos_types.jl b/src/geos_types.jl
index 5f273f2..b2506e2 100644
--- a/src/geos_types.jl
+++ b/src/geos_types.jl
@@ -1,5 +1,47 @@
 abstract type AbstractGeometry end
 
+"""
+
+    get_context(geom::AbstractGeometry)::GEOSContext
+    get_context(geometries...)::GEOSContext
+
+Return the `GEOSContext` that `geom` belongs to.
+It is also possible to pass multiple `geometries` to this function.
+In that case it is checked, that all `geometries` share the same context
+and that shared context is returned. If contexts of some geometries differ,
+an error is thrown.
+"""
+function get_context end
+function get_context(gs::AbstractVector)::GEOSContext
+    if isempty(gs)
+        get_global_context() # is this a good idea?
+    else
+        ctx = get_context(first(gs))
+        _get_context(ctx, gs)
+    end
+end
+function get_context(g1::AbstractGeometry, g2, gs...)
+    ctx = get_context(g1)
+    _get_context(ctx, g2, gs...)
+end
+function _get_context(ctx::GEOSContext, gs::AbstractVector)
+    for g in gs
+        _get_context(ctx, g)
+    end
+    ctx
+end
+function _get_context(ctx::GEOSContext, g::AbstractGeometry)
+    if ctx !== get_context(g)
+        throw(ArgumentError("Objects have distinct GEOSContext."))
+    end
+    ctx
+end
+function _get_context(ctx::GEOSContext, g, gs...)
+    _get_context(ctx, g)
+    _get_context(ctx, gs...)
+    return ctx
+end
+
 mutable struct Point <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
@@ -55,7 +97,7 @@ mutable struct MultiPoint <: AbstractGeometry
                 context),
             context)
     # create a multipoint from a list of points
-    MultiPoint(points::Vector{LibGEOS.Point}, context::GEOSContext = get_global_context()) =
+    MultiPoint(points::Vector{LibGEOS.Point}, context::GEOSContext = get_context(points)) =
         MultiPoint(
             createCollection(
                 GEOS_MULTIPOINT,
@@ -168,10 +210,10 @@ mutable struct Polygon <: AbstractGeometry
         polygon
     end
     # using 1 linear ring to form polygon with no holes - linear ring will be outer boundary of polygon
-    Polygon(ring::LinearRing, context::GEOSContext = get_global_context()) =
+    Polygon(ring::LinearRing, context::GEOSContext = get_context(ring)) =
         Polygon(ring.ptr, context)
     # using multiple linear rings to form polygon with holes - exterior linear ring will be polygon boundary and list of interior linear rings will form holes
-    Polygon(exterior::LinearRing, holes::Vector{LinearRing}, context::GEOSContext = get_global_context()) = 
+    Polygon(exterior::LinearRing, holes::Vector{LinearRing}, context::GEOSContext = get_context(exterior, holes)) = 
         Polygon(
             createPolygon(exterior.ptr,
                           GEOSGeom[ring.ptr for ring in holes],
@@ -203,7 +245,7 @@ mutable struct MultiPolygon <: AbstractGeometry
     end
 
     # create multipolygon from list of Polygon objects
-    MultiPolygon(polygons::Vector{Polygon}, context::GEOSContext = get_global_context()) =
+    MultiPolygon(polygons::Vector{Polygon}, context::GEOSContext = get_context(polygons)) =
         MultiPolygon(
             createCollection(
                 GEOS_MULTIPOLYGON,
diff --git a/src/strtree.jl b/src/strtree.jl
index 49c545b..feba255 100644
--- a/src/strtree.jl
+++ b/src/strtree.jl
@@ -37,6 +37,11 @@ function destroySTRtree(obj::STRtree)
     obj.ptr = C_NULL
 end
 
+function get_context(tree::STRtree, g::Geometry)
+    ctx = get_context(tree)
+    _get_context(ctx, g)
+end
+
 """
     query(tree::STRtree, geometry; context::GEOSContext=get_global_context())
 
@@ -50,7 +55,7 @@ Returns the objects within `tree`, whose envelope intersects the envelope of `ge
 function query(
     tree::STRtree{T},
     geometry::Geometry;
-    context::GEOSContext = get_context(tree),
+    context::GEOSContext = get_context(tree, geometry),
 ) where {T}
     matches = eltype(T)[]
 


Test patch:
diff --git a/test/test_misc.jl b/test/test_misc.jl
index 2faf7a0..34cd1a5 100644
--- a/test/test_misc.jl
+++ b/test/test_misc.jl
@@ -18,3 +18,22 @@ using LibGEOS
     end
     Point(1,2,3)
 end
+
+@testset "Context mixing" begin
+    ctx1 = LibGEOS.GEOSContext()
+    ctx2 = LibGEOS.GEOSContext()
+    p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
+    p1 = LibGEOS.Polygon(p, ctx1)
+    p2 = LibGEOS.Polygon(p, ctx2)
+
+    q = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
+    q1 = LibGEOS.Polygon(q, ctx1)
+    q2 = LibGEOS.Polygon(q, ctx2)
+    @test LibGEOS.intersects(p1, q1)
+    @test LibGEOS.intersects(p1, q1, ctx1)
+    @test LibGEOS.intersects(p2, q2)
+    @test LibGEOS.intersects(p2, q2, ctx2)
+    @test_throws ArgumentError LibGEOS.intersects(p1, q2)
+    @test_throws ArgumentError LibGEOS.intersects(p2, q1)
+end
+