Repository: JuliaGeo/LibGEOS.jl
Pull request #138 (2022-09-19 20:56:34)

Expose context within constructors and operations
I have exposed the context argument in the type constructors and within the geos operations. This addresses part of issue #91. All current tests pass so there is no regression. I also added one test per type constructor. If we think we want/need tests, let me know.


Patch:
diff --git a/src/geos_operations.jl b/src/geos_operations.jl
index b228906..625113c 100644
--- a/src/geos_operations.jl
+++ b/src/geos_operations.jl
@@ -13,54 +13,59 @@ writegeom(obj::Geometry, wktwriter::WKTWriter, context::GEOSContext = _context)
     _writegeom(obj.ptr, wktwriter, context)
 writegeom(obj::Geometry, wkbwriter::WKBWriter, context::GEOSContext = _context) =
     _writegeom(obj.ptr, wkbwriter, context)
-writegeom(obj::Geometry, context::GEOSContext = _context) = _writegeom(obj.ptr, context)
-
-function geomFromGEOS(ptr::GEOSGeom)
-    if geomTypeId(ptr) == GEOS_POINT
-        return Point(ptr)
-    elseif geomTypeId(ptr) == GEOS_LINESTRING
-        return LineString(ptr)
-    elseif geomTypeId(ptr) == GEOS_LINEARRING
-        return LinearRing(ptr)
-    elseif geomTypeId(ptr) == GEOS_POLYGON
-        return Polygon(ptr)
-    elseif geomTypeId(ptr) == GEOS_MULTIPOINT
-        return MultiPoint(ptr)
-    elseif geomTypeId(ptr) == GEOS_MULTILINESTRING
-        return MultiLineString(ptr)
-    elseif geomTypeId(ptr) == GEOS_MULTIPOLYGON
-        return MultiPolygon(ptr)
+writegeom(obj::Geometry, context::GEOSContext = _context) = 
+    _writegeom(obj.ptr, context)
+
+function geomFromGEOS(ptr::GEOSGeom, context::GEOSContext = _context)
+    id = geomTypeId(ptr, context)
+    if id == GEOS_POINT
+        return Point(ptr, context)
+    elseif id == GEOS_LINESTRING
+        return LineString(ptr, context)
+    elseif id == GEOS_LINEARRING
+        return LinearRing(ptr, context)
+    elseif id == GEOS_POLYGON
+        return Polygon(ptr, context)
+    elseif id == GEOS_MULTIPOINT
+        return MultiPoint(ptr, context)
+    elseif id == GEOS_MULTILINESTRING
+        return MultiLineString(ptr, context)
+    elseif id == GEOS_MULTIPOLYGON
+        return MultiPolygon(ptr, context)
     else
-        @assert geomTypeId(ptr) == GEOS_GEOMETRYCOLLECTION
-        return GeometryCollection(ptr)
+        @assert id == GEOS_GEOMETRYCOLLECTION
+        return GeometryCollection(ptr, context)
     end
 end
 
 readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = _context) =
-    geomFromGEOS(_readgeom(wktstring, wktreader, context))
+    geomFromGEOS(_readgeom(wktstring, wktreader, context), context)
 readgeom(wktstring::String, context::GEOSContext = _context) =
     readgeom(wktstring, WKTReader(context), context)
 
 readgeom(wkbbuffer::Vector{Cuchar}, wkbreader::WKBReader, context::GEOSContext = _context) =
-    geomFromGEOS(_readgeom(wkbbuffer, wkbreader, context))
+    geomFromGEOS(_readgeom(wkbbuffer, wkbreader, context), context)
 readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = _context) =
     readgeom(wkbbuffer, WKBReader(context), context)
 
 # -----
 # Linear referencing functions -- there are more, but these are probably sufficient for most purposes
 # -----
-project(line::LineString, point::Point) = project(line.ptr, point.ptr)
-projectNormalized(line::LineString, point::Point) = projectNormalized(line.ptr, point.ptr)
-interpolate(line::LineString, dist::Real) = Point(interpolate(line.ptr, dist))
-interpolateNormalized(line::LineString, dist::Real) =
-    Point(interpolateNormalized(line.ptr, dist))
+project(line::LineString, point::Point, context::GEOSContext = _context) =
+    project(line.ptr, point.ptr, context)
+projectNormalized(line::LineString, point::Point, context::GEOSContext = _context) =
+    projectNormalized(line.ptr, point.ptr, context)
+interpolate(line::LineString, dist::Real, context::GEOSContext = _context) =
+    Point(interpolate(line.ptr, dist, context), context)
+interpolateNormalized(line::LineString, dist::Real, context::GEOSContext = _context) =
+    Point(interpolateNormalized(line.ptr, dist, context), context)
 
 # # -----
 # # Topology operations
 # # -----
 
-buffer(obj::Geometry, dist::Real, quadsegs::Integer = 8) =
-    geomFromGEOS(buffer(obj.ptr, dist, quadsegs))
+buffer(obj::Geometry, dist::Real, quadsegs::Integer = 8, context::GEOSContext = _context) =
+    geomFromGEOS(buffer(obj.ptr, dist, quadsegs, context), context)
 bufferWithStyle(
     obj::Geometry,
     dist::Real;
@@ -68,25 +73,36 @@ bufferWithStyle(
     endCapStyle::GEOSBufCapStyles = GEOSBUF_CAP_ROUND,
     joinStyle::GEOSBufJoinStyles = GEOSBUF_JOIN_ROUND,
     mitreLimit::Real = 5.0,
-) = geomFromGEOS(
-    bufferWithStyle(obj.ptr, dist, quadsegs, endCapStyle, joinStyle, mitreLimit),
-)
-envelope(obj::Geometry) = geomFromGEOS(envelope(obj.ptr))
-envelope(obj::PreparedGeometry) = geomFromGEOS(envelope(obj.ownedby.ptr))
-minimumRotatedRectangle(obj::Geometry) = geomFromGEOS(minimumRotatedRectangle(obj.ptr))
-convexhull(obj::Geometry) = geomFromGEOS(convexhull(obj.ptr))
-boundary(obj::Geometry) = geomFromGEOS(boundary(obj.ptr))
-unaryUnion(obj::Geometry) = geomFromGEOS(unaryUnion(obj.ptr))
-pointOnSurface(obj::Geometry) = Point(pointOnSurface(obj.ptr))
-centroid(obj::Geometry) = Point(centroid(obj.ptr))
-node(obj::Geometry) = geomFromGEOS(node(obj.ptr))
-
-intersection(obj1::Geometry, obj2::Geometry) =
-    geomFromGEOS(intersection(obj1.ptr, obj2.ptr))
-difference(obj1::Geometry, obj2::Geometry) = geomFromGEOS(difference(obj1.ptr, obj2.ptr))
-symmetricDifference(obj1::Geometry, obj2::Geometry) =
-    geomFromGEOS(symmetricDifference(obj1.ptr, obj2.ptr))
-union(obj1::Geometry, obj2::Geometry) = geomFromGEOS(union(obj1.ptr, obj2.ptr))
+    context::GEOSContext = _context,) =
+    geomFromGEOS(bufferWithStyle(obj.ptr, dist, quadsegs, endCapStyle, joinStyle, mitreLimit, context), context)
+
+envelope(obj::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(envelope(obj.ptr, context), context)
+envelope(obj::PreparedGeometry, context::GEOSContext = _context) =
+    geomFromGEOS(envelope(obj.ownedby.ptr, context), context)
+minimumRotatedRectangle(obj::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(minimumRotatedRectangle(obj.ptr, context), context)
+convexhull(obj::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(convexhull(obj.ptr, context), context)
+boundary(obj::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(boundary(obj.ptr, context), context)
+unaryUnion(obj::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(unaryUnion(obj.ptr, context), context)
+pointOnSurface(obj::Geometry, context::GEOSContext = _context) =
+    Point(pointOnSurface(obj.ptr, context), context)
+centroid(obj::Geometry, context::GEOSContext = _context) =
+    Point(centroid(obj.ptr, context), context)
+node(obj::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(node(obj.ptr, context), context)
+
+intersection(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(intersection(obj1.ptr, obj2.ptr, context), context)
+difference(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(difference(obj1.ptr, obj2.ptr, context), context)
+symmetricDifference(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(symmetricDifference(obj1.ptr, obj2.ptr, context), context)
+union(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    geomFromGEOS(union(obj1.ptr, obj2.ptr, context), context)
 
 # # all arguments remain ownership of the caller (both Geometries and pointers)
 # function polygonize(geoms::Vector{GEOSGeom})
@@ -107,41 +123,53 @@ union(obj1::Geometry, obj2::Geometry) = geomFromGEOS(union(obj1.ptr, obj2.ptr))
 #     result
 # end
 
-simplify(obj::Geometry, tol::Real) = geomFromGEOS(simplify(obj.ptr, tol))
-topologyPreserveSimplify(obj::Geometry, tol::Real) =
-    geomFromGEOS(topologyPreserveSimplify(obj.ptr, tol))
-uniquePoints(obj::Geometry) = MultiPoint(uniquePoints(obj.ptr))
-delaunayTriangulationEdges(obj::Geometry, tol::Real = 0.0) =
-    MultiLineString(delaunayTriangulation(obj.ptr, tol, true))
-delaunayTriangulation(obj::Geometry, tol::Real = 0.0) =
-    GeometryCollection(delaunayTriangulation(obj.ptr, tol, false))
-constrainedDelaunayTriangulation(obj::Geometry) =
-    GeometryCollection(constrainedDelaunayTriangulation(obj.ptr))
+simplify(obj::Geometry, tol::Real, context::GEOSContext = _context) =
+    geomFromGEOS(simplify(obj.ptr, tol, context), context)
+topologyPreserveSimplify(obj::Geometry, tol::Real, context::GEOSContext = _context) =
+    geomFromGEOS(topologyPreserveSimplify(obj.ptr, tol, context), context)
+uniquePoints(obj::Geometry, context::GEOSContext = _context) =
+    MultiPoint(uniquePoints(obj.ptr, context), context)
+delaunayTriangulationEdges(obj::Geometry, tol::Real = 0.0, context::GEOSContext = _context) =
+    MultiLineString(delaunayTriangulation(obj.ptr, tol, true, context), context)
+delaunayTriangulation(obj::Geometry, tol::Real = 0.0, context::GEOSContext = _context) =
+    GeometryCollection(delaunayTriangulation(obj.ptr, tol, false, context), context)
+constrainedDelaunayTriangulation(obj::Geometry, context::GEOSContext = _context) =
+    GeometryCollection(constrainedDelaunayTriangulation(obj.ptr, context), context)
 
 
-sharedPaths(obj1::LineString, obj2::LineString) =
-    GeometryCollection(sharedPaths(obj1.ptr, obj2.ptr))
+sharedPaths(obj1::LineString, obj2::LineString, context::GEOSContext = _context) =
+    GeometryCollection(sharedPaths(obj1.ptr, obj2.ptr, context), context)
 
 # # Snap first geometry on to second with given tolerance
-snap(obj1::Geometry, obj2::Geometry, tol::Real) =
-    geomFromGEOS(snap(obj1.ptr, obj2.ptr, tol))
+snap(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = _context) =
+    geomFromGEOS(snap(obj1.ptr, obj2.ptr, tol, context), context)
 
 # -----
 # Binary predicates
 # -----
 
-disjoint(obj1::Geometry, obj2::Geometry) = disjoint(obj1.ptr, obj2.ptr)
-touches(obj1::Geometry, obj2::Geometry) = touches(obj1.ptr, obj2.ptr)
-intersects(obj1::Geometry, obj2::Geometry) = intersects(obj1.ptr, obj2.ptr)
-crosses(obj1::Geometry, obj2::Geometry) = crosses(obj1.ptr, obj2.ptr)
-within(obj1::Geometry, obj2::Geometry) = within(obj1.ptr, obj2.ptr)
-Base.contains(obj1::Geometry, obj2::Geometry) = Base.contains(obj1.ptr, obj2.ptr)
-overlaps(obj1::Geometry, obj2::Geometry) = overlaps(obj1.ptr, obj2.ptr)
-equals(obj1::Geometry, obj2::Geometry) = equals(obj1.ptr, obj2.ptr)
-equalsexact(obj1::Geometry, obj2::Geometry, tol::Real) =
-    equalsexact(obj1.ptr, obj2.ptr, tol)
-covers(obj1::Geometry, obj2::Geometry) = covers(obj1.ptr, obj2.ptr)
-coveredby(obj1::Geometry, obj2::Geometry) = coveredby(obj1.ptr, obj2.ptr)
+disjoint(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    disjoint(obj1.ptr, obj2.ptr, context)
+touches(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    touches(obj1.ptr, obj2.ptr, context)
+intersects(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    intersects(obj1.ptr, obj2.ptr, context)
+crosses(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    crosses(obj1.ptr, obj2.ptr, context)
+within(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    within(obj1.ptr, obj2.ptr, context)
+Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    Base.contains(obj1.ptr, obj2.ptr, context)
+overlaps(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    overlaps(obj1.ptr, obj2.ptr, context)
+equals(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    equals(obj1.ptr, obj2.ptr, context)
+equalsexact(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = _context) =
+    equalsexact(obj1.ptr, obj2.ptr, tol, context)
+covers(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    covers(obj1.ptr, obj2.ptr, context)
+coveredby(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    coveredby(obj1.ptr, obj2.ptr, context)
 
 
 # # -----
@@ -150,17 +178,26 @@ coveredby(obj1::Geometry, obj2::Geometry) = coveredby(obj1.ptr, obj2.ptr)
 
 prepareGeom(obj::Geometry, context::GEOSContext = _context) =
     PreparedGeometry(prepareGeom(obj.ptr, context), obj)
-Base.contains(obj1::PreparedGeometry, obj2::Geometry) = prepcontains(obj1.ptr, obj2.ptr)
-containsproperly(obj1::PreparedGeometry, obj2::Geometry) =
-    prepcontainsproperly(obj1.ptr, obj2.ptr)
-coveredby(obj1::PreparedGeometry, obj2::Geometry) = prepcoveredby(obj1.ptr, obj2.ptr)
-covers(obj1::PreparedGeometry, obj2::Geometry) = prepcovers(obj1.ptr, obj2.ptr)
-crosses(obj1::PreparedGeometry, obj2::Geometry) = prepcrosses(obj1.ptr, obj2.ptr)
-disjoint(obj1::PreparedGeometry, obj2::Geometry) = prepdisjoint(obj1.ptr, obj2.ptr)
-intersects(obj1::PreparedGeometry, obj2::Geometry) = prepintersects(obj1.ptr, obj2.ptr)
-overlaps(obj1::PreparedGeometry, obj2::Geometry) = prepoverlaps(obj1.ptr, obj2.ptr)
-touches(obj1::PreparedGeometry, obj2::Geometry) = preptouches(obj1.ptr, obj2.ptr)
-within(obj1::PreparedGeometry, obj2::Geometry) = prepwithin(obj1.ptr, obj2.ptr)
+Base.contains(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepcontains(obj1.ptr, obj2.ptr, context)
+containsproperly(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepcontainsproperly(obj1.ptr, obj2.ptr, context)
+coveredby(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepcoveredby(obj1.ptr, obj2.ptr, context)
+covers(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context)= 
+    prepcovers(obj1.ptr, obj2.ptr, context)
+crosses(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepcrosses(obj1.ptr, obj2.ptr, context)
+disjoint(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepdisjoint(obj1.ptr, obj2.ptr, context)
+intersects(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepintersects(obj1.ptr, obj2.ptr, context)
+overlaps(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepoverlaps(obj1.ptr, obj2.ptr, context)
+touches(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    preptouches(obj1.ptr, obj2.ptr, context)
+within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+    prepwithin(obj1.ptr, obj2.ptr, context)
 
 # # -----
 # # STRtree functions
@@ -175,14 +212,21 @@ within(obj1::PreparedGeometry, obj2::Geometry) = prepwithin(obj1.ptr, obj2.ptr)
 # # -----
 # # Unary predicate - return 2 on exception, 1 on true, 0 on false
 # # -----
-isEmpty(obj::Geometry) = isEmpty(obj.ptr)
-isEmpty(obj::PreparedGeometry) = isEmpty(obj.ownedby.ptr)
-isSimple(obj::Geometry) = isSimple(obj.ptr)
-isRing(obj::Geometry) = isRing(obj.ptr)
-isValid(obj::Geometry) = isValid(obj.ptr)
-hasZ(obj::Geometry) = hasZ(obj.ptr)
-
-isClosed(obj::LineString) = isClosed(obj.ptr) # Call only on LINESTRING
+isEmpty(obj::Geometry, context::GEOSContext = _context) =
+    isEmpty(obj.ptr, context)
+isEmpty(obj::PreparedGeometry, context::GEOSContext = _context) =
+    isEmpty(obj.ownedby.ptr, context)
+isSimple(obj::Geometry, context::GEOSContext = _context) =
+    isSimple(obj.ptr, context)
+isRing(obj::Geometry, context::GEOSContext = _context) =
+    isRing(obj.ptr, context)
+isValid(obj::Geometry, context::GEOSContext = _context) =
+    isValid(obj.ptr, context)
+hasZ(obj::Geometry, context::GEOSContext = _context) =
+    hasZ(obj.ptr, context)
+
+isClosed(obj::LineString, context::GEOSContext = _context) =
+    isClosed(obj.ptr, context) # Call only on LINESTRING
 
 # # -----
 # # Dimensionally Extended 9 Intersection Model related
@@ -234,16 +278,23 @@ numGeometries(obj::Geometry, context::GEOSContext = _context) =
 # 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 = _context) =
-    geomFromGEOS(getGeometry(obj.ptr, n, context))
+    geomFromGEOS(getGeometry(obj.ptr, n, context), context)
 getGeometries(obj::Geometry, context::GEOSContext = _context) =
-    geomFromGEOS.(getGeometries(obj.ptr, context))
+    [geomFromGEOS(gptr, context) for gptr in getGeometries(obj.ptr, context)]
 
 # Converts Geometry to normal form (or canonical form).
-normalize!(obj::Geometry) = normalize!(obj.ptr)
-
-interiorRing(obj::Polygon, n::Integer) = LinearRing(interiorRing(obj.ptr, n))
-interiorRings(obj::Polygon) = map(LinearRing, interiorRings(obj.ptr))
-exteriorRing(obj::Polygon) = LinearRing(exteriorRing(obj.ptr))
+normalize!(obj::Geometry, context::GEOSContext = _context) =
+    normalize!(obj.ptr, context)
+
+# LinearRings in Polygons
+numInteriorRings(obj::Polygon, context::GEOSContext = _context) =
+    numInteriorRings(obj.ptr, context)
+interiorRing(obj::Polygon, n::Integer, context::GEOSContext = _context) =
+    LinearRing(interiorRing(obj.ptr, n, context), context)
+interiorRings(obj::Polygon, context::GEOSContext = _context) =
+    map(LinearRing, interiorRings(obj.ptr, context))
+exteriorRing(obj::Polygon, context::GEOSContext = _context) =
+    LinearRing(exteriorRing(obj.ptr, context), context)
 
 # # Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
 # function getCoordSeq(ptr::GEOSGeom)
@@ -270,30 +321,39 @@ exteriorRing(obj::Polygon) = LinearRing(exteriorRing(obj.ptr))
 #     result
 # end
 
-numPoints(obj::LineString) = numPoints(obj.ptr) # Call only on LINESTRING
-startPoint(obj::LineString) = Point(startPoint(obj.ptr)) # Call only on LINESTRING
-endPoint(obj::LineString) = Point(endPoint(obj.ptr)) # Call only on LINESTRING
+numPoints(obj::LineString, context::GEOSContext = _context) =
+    numPoints(obj.ptr, context) # Call only on LINESTRING
+startPoint(obj::LineString, context::GEOSContext = _context) =
+    Point(startPoint(obj.ptr, context), context) # Call only on LINESTRING
+endPoint(obj::LineString, context::GEOSContext = _context) =
+    Point(endPoint(obj.ptr, context), context) # Call only on LINESTRING
 
 # # -----
 # # Misc functions
 # # -----
 
-area(obj::Geometry) = geomArea(obj.ptr)
-geomLength(obj::Geometry) = geomLength(obj.ptr)
+area(obj::Geometry, context::GEOSContext = _context) =
+    geomArea(obj.ptr, context)
+geomLength(obj::Geometry, context::GEOSContext = _context) =
+    geomLength(obj.ptr, context)
+
+distance(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+    geomDistance(obj1.ptr, obj2.ptr, context)
+hausdorffdistance(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) = 
+    hausdorffdistance(obj1.ptr, obj2.ptr, context)
 
-distance(obj1::Geometry, obj2::Geometry) = geomDistance(obj1.ptr, obj2.ptr)
-hausdorffdistance(obj1::Geometry, obj2::Geometry) = hausdorffdistance(obj1.ptr, obj2.ptr)
-hausdorffdistance(obj1::Geometry, obj2::Geometry, densify::Real) =
-    hausdorffdistance(obj1.ptr, obj2.ptr, densify)
+hausdorffdistance(obj1::Geometry, obj2::Geometry, densify::Real,context::GEOSContext = _context) =
+    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)
-    points = nearestPoints(obj1.ptr, obj2.ptr)
+function nearestPoints(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context)
+    points = nearestPoints(obj1.ptr, obj2.ptr, context)
     if points == C_NULL
         return Point[]
     else
-        return Point[Point(getCoordinates(points, 1)), Point(getCoordinates(points, 2))]
+        return Point[Point(getCoordinates(points, 1, context), context), 
+                     Point(getCoordinates(points, 2, context), context)]
     end
 end
 
@@ -307,8 +367,8 @@ setPrecision(
     obj::Geometry,
     grid::Real;
     flags = GEOS_PREC_VALID_OUTPUT,
-    context::GEOSContext = _context,
-) = setPrecision(obj.ptr, grid::Real, flags, context)
+    context::GEOSContext = _context,) =
+    setPrecision(obj.ptr, grid::Real, flags, context)
 
 # ----
 #  Geometry information functions
diff --git a/src/geos_types.jl b/src/geos_types.jl
index 017fc50..cd6f37b 100644
--- a/src/geos_types.jl
+++ b/src/geos_types.jl
@@ -3,70 +3,80 @@ abstract type AbstractGeometry end
 mutable struct Point <: AbstractGeometry
     ptr::GEOSGeom
     # create a point from a pointer - only makese sense if it is a pointer to a point, otherwise error
-    function Point(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function Point(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         point = if id == GEOS_POINT
-            point = new(cloneGeom(ptr))
+            point = new(cloneGeom(ptr, context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a point (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a point (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, point)
         point
     end
     # create a point from a vector of floats
-    Point(coords::Vector{Float64}) = Point(createPoint(coords))
-    Point(x::Real, y::Real) = Point(createPoint(x, y))
-    Point(x::Real, y::Real, z::Real) = Point(createPoint(x, y, z))
+    Point(coords::Vector{Float64}, context::GEOSContext = _context) =
+        Point(createPoint(coords, context), context)
+    Point(x::Real, y::Real, context::GEOSContext = _context) =
+        Point(createPoint(x, y, context), context)
+    Point(x::Real, y::Real, z::Real, context::GEOSContext = _context) =
+        Point(createPoint(x, y, z, context), context)
 end
 
 mutable struct MultiPoint <: AbstractGeometry
     ptr::GEOSGeom
-    # create a multipoint from a pointer - only makes sense if it is a pointer to a multipoint or to a point, otherwise error
-    function MultiPoint(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    # create a multipoint from a pointer - only makes sense if it is a pointer to a multipoint
+    # or to a point, otherwise error
+    function MultiPoint(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         multipoint = if id == GEOS_MULTIPOINT
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         elseif id == GEOS_POINT
             new(createCollection(GEOS_MULTIPOINT, 
-                                        GEOSGeom[cloneGeom(ptr)]))
+                                 GEOSGeom[cloneGeom(ptr, context)],
+                                 context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a multipoint (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a multipoint (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, multipoint)
         multipoint
     end
     # create a multipoint frome a vector of vector coordinates
-    MultiPoint(multipoint::Vector{Vector{Float64}}) = MultiPoint(
-        createCollection(
-            GEOS_MULTIPOINT,
-            GEOSGeom[createPoint(coords) for coords in multipoint],
-        ),
-    )
+    MultiPoint(multipoint::Vector{Vector{Float64}}, context::GEOSContext = _context) =
+        MultiPoint(
+            createCollection(
+                GEOS_MULTIPOINT,
+                GEOSGeom[createPoint(coords, context) for coords in multipoint],
+                context),
+            context)
     # create a multipoint from a list of points
-    MultiPoint(points::Vector{LibGEOS.Point}) = MultiPoint(
-        createCollection(
-            GEOS_MULTIPOINT,
-            GEOSGeom[point.ptr for point in points],
-        ),
-    )
+    MultiPoint(points::Vector{LibGEOS.Point}, context::GEOSContext = _context) =
+        MultiPoint(
+            createCollection(
+                GEOS_MULTIPOINT,
+                GEOSGeom[point.ptr for point in points],
+                context),
+            context)
 end
 
 mutable struct LineString <: AbstractGeometry
     ptr::GEOSGeom
     # create a linestring from a linestring pointer, otherwise error
-    function LineString(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function LineString(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         line = if id == GEOS_LINESTRING
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a linestring (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a linestring (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, line)
         line
     end
     #create a linestring from a list of coordiantes
-    function LineString(coords::Vector{Vector{Float64}})
-        line = new(createLineString(coords))
+    function LineString(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
+        line = new(createLineString(coords, context))
         finalizer(destroyGeom, line)
         line
     end
@@ -75,46 +85,49 @@ end
 mutable struct MultiLineString <: AbstractGeometry
     ptr::GEOSGeom
     # create a multiline string from a multilinestring or a linestring pointer, else error
-    function MultiLineString(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function MultiLineString(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         multiline = if id == GEOS_MULTILINESTRING
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         elseif id == GEOS_LINESTRING
             new(createCollection(GEOS_MULTILINESTRING, 
-                                 GEOSGeom[cloneGeom(ptr)]))
+                                 GEOSGeom[cloneGeom(ptr, context)], 
+                                 context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a multi-linestring (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a multi-linestring (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, multiline)
         multiline
     end
     # create a multilinestring from a list of linestring coordiantes
-    MultiLineString(multiline::Vector{Vector{Vector{Float64}}}) = MultiLineString(
-        createCollection(
-            GEOS_MULTILINESTRING,
-            GEOSGeom[createLineString(coords) for coords in multiline],
-        ),
-    )
-
+    MultiLineString(multiline::Vector{Vector{Vector{Float64}}},context::GEOSContext = _context) =
+        MultiLineString(
+            createCollection(
+                GEOS_MULTILINESTRING,
+                GEOSGeom[createLineString(coords, context) for coords in multiline],
+                context),
+            context)
 end
 
 mutable struct LinearRing <: AbstractGeometry
     ptr::GEOSGeom
     # create a linear ring from a linear ring pointer, otherwise error
-    function LinearRing(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function LinearRing(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         ring = if id == GEOS_LINEARRING
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a linear ring (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a linear ring (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, ring)
         ring
     end
     # create linear ring from a list of coordinates - 
     # first and last coordinates must be the same
-    function LinearRing(coords::Vector{Vector{Float64}})
-        ring = new(createLinearRing(coords))
+    function LinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
+        ring = new(createLinearRing(coords, context))
         finalizer(destroyGeom, ring)
         ring
     end
@@ -124,88 +137,106 @@ end
 mutable struct Polygon <: AbstractGeometry
     ptr::GEOSGeom
     # create polygon using GEOSGeom pointer - only makes sense if pointer points to a polygon or a linear ring to start with. 
-    function Polygon(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function Polygon(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         polygon = if id == GEOS_POLYGON
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         elseif id == GEOS_LINEARRING
-            new(cloneGeom(createPolygon(ptr)))
+            new(cloneGeom(createPolygon(ptr, context), context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a polygon (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a polygon (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, polygon)
         polygon
     end
     # using vector of coordinates in following form:
     # [[exterior], [hole1], [hole2], ...] where exterior and holeN are coordinates where the first and last point are the same
-    function Polygon(coords::Vector{Vector{Vector{Float64}}})
-        exterior = createLinearRing(coords[1])
-        interiors = GEOSGeom[createLinearRing(lr) for lr in coords[2:end]]
-        polygon = new(createPolygon(exterior, interiors))
+    function Polygon(coords::Vector{Vector{Vector{Float64}}}, context::GEOSContext = _context)
+        exterior = createLinearRing(coords[1], context)
+        interiors = GEOSGeom[createLinearRing(lr, context) for lr in coords[2:end]]
+        polygon = new(createPolygon(exterior, interiors, context))
         finalizer(destroyGeom, polygon)
         polygon
     end
     # using 1 linear ring to form polygon with no holes - linear ring will be outer boundary of polygon
-    Polygon(ring::LinearRing) = Polygon(ring.ptr)
+    Polygon(ring::LinearRing, context::GEOSContext = _context) =
+        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}) = 
-        Polygon(createPolygon(exterior.ptr, 
-                              GEOSGeom[ring.ptr for ring in holes]))
+    Polygon(exterior::LinearRing, holes::Vector{LinearRing}, context::GEOSContext = _context) = 
+        Polygon(
+            createPolygon(exterior.ptr,
+                          GEOSGeom[ring.ptr for ring in holes],
+                          context),
+            context)
 end
 
 mutable struct MultiPolygon <: AbstractGeometry
     ptr::GEOSGeom
     # create multipolygon using a multipolygon or polygon pointer, else error
-    function MultiPolygon(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function MultiPolygon(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         multipolygon = if id == GEOS_MULTIPOLYGON
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         elseif id == GEOS_POLYGON
-            new(createCollection(GEOS_MULTIPOLYGON,
-                                GEOSGeom[cloneGeom(ptr)]))
+            new(createCollection(
+                    GEOS_MULTIPOLYGON,
+                    GEOSGeom[cloneGeom(ptr, context)],
+                    context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a multi-polygon (yet). Please open an issue if you think this conversion makes sense.")    
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a multi-polygon (yet).
+                   Please open an issue if you think this conversion makes sense.")    
         end
         finalizer(destroyGeom, multipolygon)
         multipolygon
     end
 
     # create multipolygon from list of Polygon objects
-    MultiPolygon(polygons::Vector{Polygon}) = MultiPolygon(
-        createCollection(
-            GEOS_MULTIPOLYGON,
-            GEOSGeom[poly.ptr for poly in polygons],))
+    MultiPolygon(polygons::Vector{Polygon}, context::GEOSContext = _context) =
+        MultiPolygon(
+            createCollection(
+                GEOS_MULTIPOLYGON,
+                GEOSGeom[poly.ptr for poly in polygons],
+                context),
+            context)
 
     # create multipolygon using list of polygon coordinates - note that each polygon can have holes as explained above in Polygon comments
-    MultiPolygon(multipolygon::Vector{Vector{Vector{Vector{Float64}}}}) = MultiPolygon(
-        createCollection(
-            GEOS_MULTIPOLYGON,
-            GEOSGeom[
-                createPolygon(
-                    createLinearRing(coords[1]),
-                    GEOSGeom[createLinearRing(c) for c in coords[2:end]],
-                ) for coords in multipolygon
-            ],
-        ),
-    )
+    MultiPolygon(multipolygon::Vector{Vector{Vector{Vector{Float64}}}}, context::GEOSContext = _context) =
+        MultiPolygon(
+            createCollection(
+                GEOS_MULTIPOLYGON,
+                GEOSGeom[
+                    createPolygon(
+                        createLinearRing(coords[1], context),
+                        GEOSGeom[createLinearRing(c, context) for c in coords[2:end]],
+                        context)
+                    for coords in multipolygon],
+                context),
+            context)
 end
 
 mutable struct GeometryCollection <: AbstractGeometry
     ptr::GEOSGeom
     # create a geometric collection from a pointer to a geometric collection, else error
-    function GeometryCollection(ptr::GEOSGeom)
-        id = LibGEOS.geomTypeId(ptr)
+    function GeometryCollection(ptr::GEOSGeom, context::GEOSContext = _context)
+        id = LibGEOS.geomTypeId(ptr, context)
         geometrycollection = if id == GEOS_GEOMETRYCOLLECTION
-            new(cloneGeom(ptr))
+            new(cloneGeom(ptr, context))
         else
-            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a geometry collection (yet). Please open an issue if you think this conversion makes sense.")
+            error("LibGEOS: Can't convert a pointer to an element with a GeomType ID of $id to a geometry collection (yet).
+                   Please open an issue if you think this conversion makes sense.")
         end
         finalizer(destroyGeom, geometrycollection)
         geometrycollection
     end
     # create a geometric collection from a list of pointers to geometric objects
-    GeometryCollection(collection::Vector{GEOSGeom}) =
-        GeometryCollection(createCollection(GEOS_GEOMETRYCOLLECTION, collection))
+    GeometryCollection(collection::Vector{GEOSGeom}, context::GEOSContext = _context) =
+        GeometryCollection(
+            createCollection(
+                GEOS_GEOMETRYCOLLECTION,
+                collection,
+                context),
+            context)
 end
 
 const Geometry = Union{


Test patch:
diff --git a/test/test_geos_types.jl b/test/test_geos_types.jl
index 4911655..587663b 100644
--- a/test/test_geos_types.jl
+++ b/test/test_geos_types.jl
@@ -48,10 +48,16 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         @test_throws ErrorException LibGEOS.Point(LibGEOS.Polygon(
             [[[-2.0, -2.0], [2.0, 2.0],[-2.0,2.0], [-2.0, -2.0]]]).ptr)
 
+        # Test point made with local context
+        ctx = LibGEOS.GEOSContext()
+        point_ctx = LibGEOS.Point(point2D.ptr, ctx)
+        @test LibGEOS.equals(point_ctx, point2D) 
+
         LibGEOS.destroyGeom(point2D)
         LibGEOS.destroyGeom(point3D)
         LibGEOS.destroyGeom(point_coord)
         LibGEOS.destroyGeom(point_ptr)
+        LibGEOS.destroyGeom(point_ctx)
     end
 
     @testset "MultiPoint" begin
@@ -95,6 +101,11 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         @test LibGEOS.numGeometries(point_coord1) == 1
         @test LibGEOS.equals(point_coord2, mpoint_coord)
 
+        # Test MultiPoint made with local context
+        ctx = LibGEOS.GEOSContext()
+        mpoint_ctx = LibGEOS.MultiPoint(mpoint_coord.ptr, ctx)
+        @test LibGEOS.equals(mpoint_ctx, mpoint_coord) 
+
         LibGEOS.destroyGeom(mpoint_coord)
         LibGEOS.destroyGeom(point1)
         LibGEOS.destroyGeom(point2)
@@ -102,6 +113,7 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         LibGEOS.destroyGeom(mpoint_ptr)
         LibGEOS.destroyGeom(point_coord1)
         LibGEOS.destroyGeom(point_coord2)
+        LibGEOS.destroyGeom(mpoint_ctx)
     end
 
     @testset "LineString" begin
@@ -119,8 +131,14 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         @test_throws ErrorException LibGEOS.LineString(LibGEOS.Polygon(
             [[[-2.0, -2.0], [2.0, 2.0],[-2.0,2.0], [-2.0, -2.0]]]).ptr)
 
+        # Test LineString made with local context
+        ctx = LibGEOS.GEOSContext()
+        linestring_ctx = LibGEOS.LineString(ls_coord.ptr, ctx)
+        @test LibGEOS.equals(linestring_ctx, ls_coord) 
+
         LibGEOS.destroyGeom(ls_coord)
         LibGEOS.destroyGeom(ls_ptr)
+        LibGEOS.destroyGeom(linestring_ctx)
     end
 
     @testset "MultiLineString" begin
@@ -138,8 +156,14 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         @test_throws ErrorException LibGEOS.MultiLineString(LibGEOS.Polygon(
             [[[-2.0, -2.0], [2.0, 2.0],[-2.0,2.0], [-2.0, -2.0]]]).ptr)
 
+        # Test MultiLineString made with local context
+        ctx = LibGEOS.GEOSContext()
+        multilinestring_ctx = LibGEOS.MultiLineString(mls_coord.ptr, ctx)
+        @test LibGEOS.equals(multilinestring_ctx, mls_coord) 
+
         LibGEOS.destroyGeom(mls_coord)
         LibGEOS.destroyGeom(mls_ptr)
+        LibGEOS.destroyGeom(multilinestring_ctx)
     end
 
     @testset "LinearRing" begin
@@ -158,8 +182,14 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         @test_throws ErrorException LibGEOS.LinearRing(LibGEOS.Polygon(
             [[[-2.0, -2.0], [2.0, 2.0],[-2.0,2.0], [-2.0, -2.0]]]).ptr)
 
+        # Test LinearRing made with local context
+        ctx = LibGEOS.GEOSContext()
+        linearring_ctx = LibGEOS.LinearRing(lr_coord.ptr, ctx)
+        @test LibGEOS.equals(linearring_ctx, lr_coord)
+
         LibGEOS.destroyGeom(lr_coord)
         LibGEOS.destroyGeom(lr_ptr)
+        LibGEOS.destroyGeom(linearring_ctx)
     end
 
 
@@ -231,6 +261,11 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
              [[0.0, 0.0], [0.0, -1.0], [-1.0, -1.0],[-1.0, 0.0], [0.0, 0.0]]]
         @test length(LibGEOS.interiorRings(poly_rings2)) == 2
 
+        # Test Polygon made with local context
+        ctx = LibGEOS.GEOSContext()
+        poly_ctx = LibGEOS.Polygon(poly_vec.ptr, ctx)
+        @test LibGEOS.equals(poly_ctx, poly_vec)
+
         LibGEOS.destroyGeom(poly_vec)
         LibGEOS.destroyGeom(poly_ptr)
         LibGEOS.destroyGeom(ring_ext)
@@ -241,6 +276,7 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         LibGEOS.destroyGeom(ring_int2)
         LibGEOS.destroyGeom(poly_rings1)
         LibGEOS.destroyGeom(poly_rings2)
+        LibGEOS.destroyGeom(poly_ctx)
     end
 
     @testset "MultiPolygon" begin
@@ -274,12 +310,18 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         testValidTypeDims(mpoly_polylist)
         @test LibGEOS.equals(mpoly_polylist, mpoly_coord)
 
+        # Test MultiPolygon made with local context
+        ctx = LibGEOS.GEOSContext()
+        mpoly_ctx = LibGEOS.MultiPolygon(mpoly_coord.ptr, ctx)
+        @test LibGEOS.equals(mpoly_ctx, mpoly_coord)
+
         LibGEOS.destroyGeom(poly1)
         LibGEOS.destroyGeom(poly2)
         LibGEOS.destroyGeom(mpoly_coord)
         LibGEOS.destroyGeom(mpoly_ptr)
         LibGEOS.destroyGeom(mpoly_poly_ptr)
         LibGEOS.destroyGeom(mpoly_polylist)
+        LibGEOS.destroyGeom(mpoly_ctx)
     end
 
     @testset "GeometryCollections" begin
@@ -296,10 +338,16 @@ testValidTypeDims(multipoly::LibGEOS.MultiPolygon) =
         @test LibGEOS.geomTypeId(geomcollection_ptr.ptr) == LibGEOS.GEOS_GEOMETRYCOLLECTION
         @test LibGEOS.equals(geomcol_ptr_list, geomcollection_ptr)
 
+        # Test GeomertyCollections made with local context
+        ctx = LibGEOS.GEOSContext()
+        geomcollect_ctx = LibGEOS.GeometryCollection(geomcol_ptr_list.ptr, ctx)
+        @test LibGEOS.equals(geomcollect_ctx, geomcol_ptr_list)
+
         LibGEOS.destroyGeom(point)
         LibGEOS.destroyGeom(poly)
         LibGEOS.destroyGeom(geomcol_ptr_list)
         LibGEOS.destroyGeom(geomcollection_ptr)
+        LibGEOS.destroyGeom(geomcollect_ctx)
     end
 
 end
\ No newline at end of file