Repository: JuliaGeo/LibGEOS.jl
Pull request #141 (2022-10-24 21:01:35)

control global context
When writing multi threaded code, it is important, to never use the global context. To make this easier to test and debug, this PR adds the following helper function:
```julia
using LibGEOS

LibGEOS.allow_global_context!(false) do
   Point(1,2,3) # oups implicit global context
end
```
```
ERROR: LibGEOS global context disallowed, a GEOSContext must be passed explicitly.
```


Patch:
diff --git a/src/LibGEOS.jl b/src/LibGEOS.jl
index 7e27d4a..9040570 100644
--- a/src/LibGEOS.jl
+++ b/src/LibGEOS.jl
@@ -107,6 +107,20 @@ function geosjl_errorhandler(message::Ptr{UInt8}, userdata)
     end
 end
 
+"""
+
+    GEOSContext
+
+Every LibGEOS object needs to live somewhere in memory. Also many LibGEOS functions
+need scratch memory or caches to do their job.
+
+A `GEOSContext` governs such memory. Almost every function in LibGEOS accepts a `context`
+argument, that allows passing a context explicitly. If no context is passed, a global context is used.
+
+Using the global context is fine, as long as no multi threading is used. 
+If multi threading is used, the global context should be avoided and every operation should only 
+involve objects that live in the context passed to the operation.
+"""
 mutable struct GEOSContext
     ptr::Ptr{Cvoid}  # GEOSContextHandle_t
 
@@ -123,7 +137,7 @@ mutable struct GEOSContext
 end
 
 "Get a copy of a string from GEOS, freeing the GEOS managed memory."
-function string_copy_free(s::Cstring, context::Ptr{Cvoid} = _context.ptr)::String
+function string_copy_free(s::Cstring, context::Ptr{Cvoid} = get_global_context().ptr)::String
     copy = unsafe_string(s)
     GEOSFree_r(context, pointer(s))
     return copy
@@ -191,8 +205,50 @@ mutable struct WKBWriter
     end
 end
 
+const _GLOBAL_CONTEXT = Ref{GEOSContext}()
+const _GLOBAL_CONTEXT_ALLOWED = Ref(false)
+
+function get_global_context()::GEOSContext
+    if _GLOBAL_CONTEXT_ALLOWED[]
+        _GLOBAL_CONTEXT[]
+    else
+        msg = """
+        LibGEOS global context disallowed, a `GEOSContext` must be passed explicitly.
+        Alternatively you can allow the global context by calling:
+        `LibGEOS.allow_global_context!(true)`
+        """
+        error(msg)
+    end
+end
+
+"""
+
+    allow_global_context!(bool::Bool)
+
+Allow (bool=true) or disallow (bool=false) using the global LibGEOS context.
+
+
+    allow_global_context!(f, bool::Bool)
+
+Call `f` with global context usage allowed according to `bool`
+
+Generally this function should only be used as a debugging tool, mostly for multithreaded programs.
+See also [`GEOSContext`](@ref).
+"""
+function allow_global_context!(bool::Bool)
+    _GLOBAL_CONTEXT_ALLOWED[] = bool
+end
+
+function allow_global_context!(f, bool::Bool)
+    old = _GLOBAL_CONTEXT_ALLOWED[]
+    allow_global_context!(bool)
+    f()
+    allow_global_context!(old)
+end
+
 function __init__()
-    global _context = GEOSContext()
+    _GLOBAL_CONTEXT_ALLOWED[] = true
+    _GLOBAL_CONTEXT[] = GEOSContext()
 end
 
 include("geos_functions.jl")
diff --git a/src/geos_functions.jl b/src/geos_functions.jl
index dfffbad..b4e627a 100644
--- a/src/geos_functions.jl
+++ b/src/geos_functions.jl
@@ -1,17 +1,17 @@
-function _readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = _context)
+function _readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = get_global_context())
     result = GEOSWKTReader_read_r(context.ptr, wktreader.ptr, pointer(wktstring))
     if result == C_NULL
         error("LibGEOS: Error in GEOSWKTReader_read_r while reading $wktstring")
     end
     result
 end
-_readgeom(wktstring::String, context::GEOSContext = _context) =
+_readgeom(wktstring::String, context::GEOSContext = get_global_context()) =
     _readgeom(wktstring, WKTReader(context), context)
 
 function _readgeom(
     wkbbuffer::Vector{Cuchar},
     wkbreader::WKBReader,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSWKBReader_read_r(
         context.ptr,
@@ -24,20 +24,20 @@ function _readgeom(
     end
     result
 end
-_readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = _context) =
+_readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = get_global_context()) =
     _readgeom(wkbbuffer, WKBReader(context), context)
 
-function _writegeom(geom::GEOSGeom, wktwriter::WKTWriter, context::GEOSContext = _context)
+function _writegeom(geom::GEOSGeom, wktwriter::WKTWriter, context::GEOSContext = get_global_context())
     GEOSWKTWriter_write_r(context.ptr, wktwriter.ptr, geom)
 end
 
-function _writegeom(geom::GEOSGeom, wkbwriter::WKBWriter, context::GEOSContext = _context)
+function _writegeom(geom::GEOSGeom, wkbwriter::WKBWriter, context::GEOSContext = get_global_context())
     wkbsize = Ref{Csize_t}()
     result = GEOSWKBWriter_write_r(context.ptr, wkbwriter.ptr, geom, wkbsize)
     unsafe_wrap(Array, result, wkbsize[], own = true)
 end
 
-_writegeom(geom::GEOSGeom, context::GEOSContext = _context) =
+_writegeom(geom::GEOSGeom, context::GEOSContext = get_global_context()) =
     _writegeom(geom, WKTWriter(context), context)
 
 # -----
@@ -48,7 +48,7 @@ _writegeom(geom::GEOSGeom, context::GEOSContext = _context) =
 
 Create a Coordinate sequence with ``size'' coordinates of ``dims'' dimensions (Return NULL on exception)
 """
-function createCoordSeq(size::Integer, context::GEOSContext = _context; ndim::Integer = 2)
+function createCoordSeq(size::Integer, context::GEOSContext = get_global_context(); ndim::Integer = 2)
     @assert ndim >= 2
     result = GEOSCoordSeq_create_r(context.ptr, size, ndim)
     if result == C_NULL
@@ -58,7 +58,7 @@ function createCoordSeq(size::Integer, context::GEOSContext = _context; ndim::In
 end
 
 # Clone a Coordinate Sequence (Return NULL on exception)
-function cloneCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function cloneCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     result = GEOSCoordSeq_clone_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSCoordSeq_clone")
@@ -66,7 +66,7 @@ function cloneCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
     result
 end
 
-function destroyCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function destroyCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     result = GEOSCoordSeq_destroy_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSCoordSeq_destroy")
@@ -75,7 +75,7 @@ function destroyCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
 end
 
 # Set ordinate values in a Coordinate Sequence (Return 0 on exception)
-function setX!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
+function setX!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -88,7 +88,7 @@ function setX!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext
     result
 end
 
-function setY!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
+function setY!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -101,7 +101,7 @@ function setY!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext
     result
 end
 
-function setZ!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
+function setZ!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -115,7 +115,7 @@ function setZ!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext
 end
 
 "Get size info from a Coordinate Sequence (Return 0 on exception)"
-function getSize(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function getSize(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     out = Ref{UInt32}()
     result = GEOSCoordSeq_getSize_r(context.ptr, ptr, out)
     if result == 0
@@ -125,7 +125,7 @@ function getSize(ptr::GEOSCoordSeq, context::GEOSContext = _context)
 end
 
 "Get dimensions info from a Coordinate Sequence (Return 0 on exception)"
-function getDimensions(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function getDimensions(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     out = Ref{UInt32}()
     result = GEOSCoordSeq_getDimensions_r(context.ptr, ptr, out)
     if result == 0
@@ -139,7 +139,7 @@ function setCoordSeq!(
     ptr::GEOSCoordSeq,
     i::Integer,
     coords::Vector{Float64},
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     ndim = length(coords)
     @assert ndim >= 2
@@ -159,7 +159,7 @@ end
 
 Create a createCoordSeq of a single 2D coordinate
 """
-function createCoordSeq(x::Real, y::Real, context::GEOSContext = _context)
+function createCoordSeq(x::Real, y::Real, context::GEOSContext = get_global_context())
     coordinates = createCoordSeq(1, context, ndim = 2)
     setX!(coordinates, 1, x, context)
     setY!(coordinates, 1, y, context)
@@ -171,7 +171,7 @@ end
 
 Create a createCoordSeq of a single 3D coordinate
 """
-function createCoordSeq(x::Real, y::Real, z::Real, context::GEOSContext = _context)
+function createCoordSeq(x::Real, y::Real, z::Real, context::GEOSContext = get_global_context())
     coordinates = createCoordSeq(1, context, ndim = 3)
     setX!(coordinates, 1, x, context)
     setY!(coordinates, 1, y, context)
@@ -184,7 +184,7 @@ end
 
 Create a createCoordSeq of a single N dimensional coordinate
 """
-function createCoordSeq(coords::Vector{Float64}, context::GEOSContext = _context)
+function createCoordSeq(coords::Vector{Float64}, context::GEOSContext = get_global_context())
     ndim = length(coords)
     @assert ndim >= 2
     coordinates = createCoordSeq(1, context, ndim = ndim)
@@ -196,7 +196,7 @@ end
 
 Create a createCoordSeq of a multiple N dimensional coordinate
 """
-function createCoordSeq(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
+function createCoordSeq(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context())
     ncoords = length(coords)
     @assert ncoords > 0
     ndim = length(coords[1])
@@ -207,7 +207,7 @@ function createCoordSeq(coords::Vector{Vector{Float64}}, context::GEOSContext =
     coordinates
 end
 
-function getX(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
+function getX(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -221,7 +221,7 @@ function getX(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
     out[]
 end
 
-function getX(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function getX(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     ncoords = getSize(ptr, context)
     xcoords = Array{Float64}(ncoords)
     start = pointer(xcoords)
@@ -232,7 +232,7 @@ function getX(ptr::GEOSCoordSeq, context::GEOSContext = _context)
     xcoords
 end
 
-function getY(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
+function getY(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -246,7 +246,7 @@ function getY(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
     out[]
 end
 
-function getY(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function getY(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     ncoords = getSize(ptr, context)
     ycoords = Array{Float64}(ncoords)
     start = pointer(ycoords)
@@ -257,7 +257,7 @@ function getY(ptr::GEOSCoordSeq, context::GEOSContext = _context)
     ycoords
 end
 
-function getZ(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
+function getZ(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -271,7 +271,7 @@ function getZ(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
     out[]
 end
 
-function getZ(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function getZ(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     ncoords = getSize(ptr, context)
     zcoords = Array{Float64}(ncoords)
     start = pointer(zcoords)
@@ -282,7 +282,7 @@ function getZ(ptr::GEOSCoordSeq, context::GEOSContext = _context)
     zcoords
 end
 
-function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
+function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = get_global_context())
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
@@ -300,7 +300,7 @@ function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _c
     coord
 end
 
-function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     ndim = getDimensions(ptr, context)
     ncoords = getSize(ptr, context)
     coordseq = Vector{Float64}[]
@@ -311,7 +311,7 @@ function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = _context)
     coordseq
 end
 
-function isCCW(ptr::GEOSCoordSeq, context::GEOSContext = _context)::Bool
+function isCCW(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())::Bool
     d = UInt8[1]
     GC.@preserve result = GEOSCoordSeq_isCCW_r(context.ptr, ptr, pointer(d))
     if result == C_NULL
@@ -326,11 +326,11 @@ end
 # (GEOSGeometry ownership is retained by caller)
 
 # Return distance of point 'p' projected on 'g' from origin of 'g'. Geometry 'g' must be a lineal geometry
-project(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = _context) =
+project(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = get_global_context()) =
     GEOSProject_r(context.ptr, g, p)
 
 # Return closest point to given distance within geometry (Geometry must be a LineString)
-function interpolate(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
+function interpolate(ptr::GEOSGeom, d::Real, context::GEOSContext = get_global_context())
     result = GEOSInterpolate_r(context.ptr, ptr, d)
     if result == C_NULL
         error("LibGEOS: Error in GEOSInterpolate")
@@ -338,10 +338,10 @@ function interpolate(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
     result
 end
 
-projectNormalized(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = _context) =
+projectNormalized(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = get_global_context()) =
     GEOSProjectNormalized_r(context.ptr, g, p)
 
-function interpolateNormalized(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
+function interpolateNormalized(ptr::GEOSGeom, d::Real, context::GEOSContext = get_global_context())
     result = GEOSInterpolateNormalized_r(context.ptr, ptr, d)
     if result == C_NULL
         error("LibGEOS: Error in GEOSInterpolateNormalized")
@@ -357,7 +357,7 @@ end
 # The user can control the accuracy of the curve approximation by specifying the number of linear segments with which to approximate a curve.
 
 # Always returns a polygon. The negative or zero-distance buffer of lines and points is always an empty Polygon.
-buffer(ptr::GEOSGeom, width::Real, quadsegs::Integer = 8, context::GEOSContext = _context) =
+buffer(ptr::GEOSGeom, width::Real, quadsegs::Integer = 8, context::GEOSContext = get_global_context()) =
     GEOSBuffer_r(context.ptr, ptr, width, Int32(quadsegs))
 
 bufferWithStyle(
@@ -367,7 +367,7 @@ bufferWithStyle(
     endCapStyle::GEOSBufCapStyles = GEOSBUF_CAP_ROUND,
     joinStyle::GEOSBufJoinStyles = GEOSBUF_JOIN_ROUND,
     mitreLimit::Real = 5.0,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 ) = GEOSBufferWithStyle_r(
     context.ptr,
     ptr,
@@ -393,40 +393,40 @@ bufferWithStyle(
 # -----
 # (GEOSCoordSequence* arguments will become ownership of the returned object.)
 
-function createPoint(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function createPoint(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     result = GEOSGeom_createPoint_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_createPoint")
     end
     result
 end
-createPoint(x::Real, y::Real, context::GEOSContext = _context) =
+createPoint(x::Real, y::Real, context::GEOSContext = get_global_context()) =
     createPoint(createCoordSeq(x, y, context), context)
-createPoint(x::Real, y::Real, z::Real, context::GEOSContext = _context) =
+createPoint(x::Real, y::Real, z::Real, context::GEOSContext = get_global_context()) =
     createPoint(createCoordSeq(x, y, z, context), context)
-createPoint(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
+createPoint(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context()) =
     createPoint(createCoordSeq(coords, context), context)
-createPoint(coords::Vector{Float64}, context::GEOSContext = _context) =
+createPoint(coords::Vector{Float64}, context::GEOSContext = get_global_context()) =
     createPoint(createCoordSeq(Vector{Float64}[coords], context), context)
 
-function createLinearRing(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function createLinearRing(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     result = GEOSGeom_createLinearRing_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_createLinearRing")
     end
     result
 end
-createLinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
+createLinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context()) =
     GEOSGeom_createLinearRing_r(context.ptr, createCoordSeq(coords, context))
 
-function createLineString(ptr::GEOSCoordSeq, context::GEOSContext = _context)
+function createLineString(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     result = GEOSGeom_createLineString_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_createLineString")
     end
     result
 end
-createLineString(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
+createLineString(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context()) =
     GEOSGeom_createLineString_r(context.ptr, createCoordSeq(coords, context))
 
 # Second argument is an array of GEOSGeometry* objects.
@@ -435,7 +435,7 @@ createLineString(coords::Vector{Vector{Float64}}, context::GEOSContext = _contex
 function createPolygon(
     shell::GEOSGeom,
     holes::Vector{GEOSGeom},
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSGeom_createPolygon_r(context.ptr, shell, pointer(holes), length(holes))
     if result == C_NULL
@@ -444,13 +444,13 @@ function createPolygon(
     result
 end
 # convenience function to create polygon without holes
-createPolygon(shell::GEOSGeom, context::GEOSContext = _context) =
+createPolygon(shell::GEOSGeom, context::GEOSContext = get_global_context()) =
     createPolygon(shell, GEOSGeom[], context)
 
 function createCollection(
     geomtype::GEOSGeomTypes,
     geoms::Vector{GEOSGeom},
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result =
         GEOSGeom_createCollection_r(context.ptr, geomtype, pointer(geoms), length(geoms))
@@ -460,7 +460,7 @@ function createCollection(
     result
 end
 
-function createEmptyCollection(geomtype::GEOSGeomTypes, context::GEOSContext = _context)
+function createEmptyCollection(geomtype::GEOSGeomTypes, context::GEOSContext = get_global_context())
     result = GEOSGeom_createEmptyCollection_r(context.ptr, geomtype)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_createEmptyCollection")
@@ -469,7 +469,7 @@ function createEmptyCollection(geomtype::GEOSGeomTypes, context::GEOSContext = _
 end
 
 # Memory management
-function cloneGeom(ptr::GEOSGeom, context::GEOSContext = _context)
+function cloneGeom(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeom_clone_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_clone")
@@ -477,7 +477,7 @@ function cloneGeom(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function destroyGeom(ptr::GEOSGeom, context::GEOSContext = _context)
+function destroyGeom(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     GC.@preserve context begin
         GEOSGeom_destroy_r(context.ptr, ptr)
     end
@@ -487,7 +487,7 @@ end
 # Topology operations - return NULL on exception.
 # -----
 
-function envelope(ptr::GEOSGeom, context::GEOSContext = _context)
+function envelope(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSEnvelope_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSEnvelope")
@@ -504,7 +504,7 @@ has width equal to the minimum diameter, and a longer length. If the convex hill
 degenerate (a line or point) a LINESTRING or POINT is returned. The minimum rotated rectangle can
 be used as an extremely generalized representation for the given geometry.
 """
-function minimumRotatedRectangle(ptr::GEOSGeom, context::GEOSContext = _context)
+function minimumRotatedRectangle(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSMinimumRotatedRectangle_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSMinimumRotatedRectangle")
@@ -513,7 +513,7 @@ function minimumRotatedRectangle(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 
-function intersection(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function intersection(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSIntersection_r(context.ptr, g1, g2)
     if result == C_NULL
         error("LibGEOS: Error in GEOSIntersection")
@@ -525,7 +525,7 @@ end
 
 # Returns:
 # if the convex hull contains 3 or more points, a Polygon; 2 points, a LineString; 1 point, a Point; 0 points, an empty GeometryCollection.
-function convexhull(ptr::GEOSGeom, context::GEOSContext = _context)
+function convexhull(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSConvexHull_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSConvexHull")
@@ -533,7 +533,7 @@ function convexhull(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function difference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function difference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSDifference_r(context.ptr, g1, g2)
     if result == C_NULL
         error("LibGEOS: Error in GEOSDifference")
@@ -541,7 +541,7 @@ function difference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function symmetricDifference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function symmetricDifference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSSymDifference_r(context.ptr, g1, g2)
     if result == C_NULL
         error("LibGEOS: Error in GEOSSymDifference")
@@ -549,7 +549,7 @@ function symmetricDifference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext =
     result
 end
 
-function boundary(ptr::GEOSGeom, context::GEOSContext = _context)
+function boundary(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSBoundary_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSBoundary")
@@ -557,7 +557,7 @@ function boundary(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function union(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function union(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSUnion_r(context.ptr, g1, g2)
     if result == C_NULL
         error("LibGEOS: Error in GEOSUnion")
@@ -565,7 +565,7 @@ function union(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function unaryUnion(ptr::GEOSGeom, context::GEOSContext = _context)
+function unaryUnion(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSUnaryUnion_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSUnaryUnion")
@@ -573,7 +573,7 @@ function unaryUnion(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function pointOnSurface(ptr::GEOSGeom, context::GEOSContext = _context)
+function pointOnSurface(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSPointOnSurface_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSPointOnSurface")
@@ -581,7 +581,7 @@ function pointOnSurface(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function centroid(ptr::GEOSGeom, context::GEOSContext = _context)
+function centroid(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGetCentroid_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGetCentroid")
@@ -589,7 +589,7 @@ function centroid(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function node(ptr::GEOSGeom, context::GEOSContext = _context)
+function node(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSNode_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSNode")
@@ -598,7 +598,7 @@ function node(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # all arguments remain ownership of the caller (both Geometries and pointers)
-function polygonize(geoms::Vector{GEOSGeom}, context::GEOSContext = _context)
+function polygonize(geoms::Vector{GEOSGeom}, context::GEOSContext = get_global_context())
     result = GEOSPolygonize_r(context.ptr, pointer(geoms), length(geoms))
     if result == C_NULL
         error("LibGEOS: Error in GEOSPolygonize")
@@ -608,7 +608,7 @@ end
 # GEOSPolygonizer_getCutEdges
 # GEOSPolygonize_full
 
-function lineMerge(ptr::GEOSGeom, context::GEOSContext = _context)
+function lineMerge(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSLineMerge_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSLineMerge")
@@ -616,7 +616,7 @@ function lineMerge(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function simplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
+function simplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = get_global_context())
     result = GEOSSimplify_r(context.ptr, ptr, tol)
     if result == C_NULL
         error("LibGEOS: Error in GEOSSimplify")
@@ -624,7 +624,7 @@ function simplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
     result
 end
 
-function topologyPreserveSimplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
+function topologyPreserveSimplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = get_global_context())
     result = GEOSTopologyPreserveSimplify_r(context.ptr, ptr, tol)
     if result == C_NULL
         error("LibGEOS: Error in GEOSTopologyPreserveSimplify")
@@ -634,7 +634,7 @@ end
 
 # Return all distinct vertices of input geometry as a MULTIPOINT.
 # (Note that only 2 dimensions of the vertices are considered when testing for equality)
-function uniquePoints(ptr::GEOSGeom, context::GEOSContext = _context)
+function uniquePoints(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeom_extractUniquePoints_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_extractUniquePoints")
@@ -649,7 +649,7 @@ end
 #  - second element is a MULTILINESTRING containing shared paths
 #    having the _opposite_ direction on the two inputs
 # (Returns NULL on exception)
-function sharedPaths(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function sharedPaths(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSSharedPaths_r(context.ptr, g1, g2)
     if result == C_NULL
         error("LibGEOS: Error in GEOSSharedPaths")
@@ -659,7 +659,7 @@ end
 
 # Snap first geometry on to second with given tolerance
 # (Returns a newly allocated geometry, or NULL on exception)
-function snap(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = _context)
+function snap(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = get_global_context())
     result = GEOSSnap_r(context.ptr, g1, g2, tol)
     if result == C_NULL
         error("LibGEOS: Error in GEOSSnap")
@@ -679,7 +679,7 @@ function delaunayTriangulation(
     ptr::GEOSGeom,
     tol::Real = 0.0,
     onlyEdges::Bool = false,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSDelaunayTriangulation_r(context.ptr, ptr, tol, Int32(onlyEdges))
     if result == C_NULL
@@ -688,7 +688,7 @@ function delaunayTriangulation(
     result
 end
 
-function constrainedDelaunayTriangulation(ptr::GEOSGeom, context::GEOSContext = _context)
+function constrainedDelaunayTriangulation(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSConstrainedDelaunayTriangulation_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSConstrainedDelaunayTriangulation")
@@ -699,7 +699,7 @@ end
 # -----
 # Binary predicates - return 2 on exception, 1 on true, 0 on false
 # -----
-function disjoint(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function disjoint(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSDisjoint_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSDisjoint")
@@ -707,7 +707,7 @@ function disjoint(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function touches(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function touches(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSTouches_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSTouches")
@@ -715,7 +715,7 @@ function touches(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function intersects(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function intersects(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSIntersects_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSIntersects")
@@ -723,7 +723,7 @@ function intersects(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function crosses(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function crosses(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSCrosses_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSCrosses")
@@ -731,7 +731,7 @@ function crosses(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function within(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function within(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSWithin_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSWithin")
@@ -739,7 +739,7 @@ function within(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function Base.contains(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function Base.contains(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSContains_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSContains")
@@ -747,7 +747,7 @@ function Base.contains(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _conte
     result != 0x00
 end
 
-function overlaps(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function overlaps(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSOverlaps_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSOverlaps")
@@ -755,7 +755,7 @@ function overlaps(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function equals(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function equals(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSEquals_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSEquals")
@@ -763,7 +763,7 @@ function equals(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function equalsexact(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = _context)
+function equalsexact(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = get_global_context())
     result = GEOSEqualsExact_r(context.ptr, g1, g2, tol)
     if result == 0x02
         error("LibGEOS: Error in GEOSEqualsExact")
@@ -771,7 +771,7 @@ function equalsexact(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext
     result != 0x00
 end
 
-function covers(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function covers(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSCovers_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSCovers")
@@ -779,7 +779,7 @@ function covers(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function coveredby(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function coveredby(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSCoveredBy_r(context.ptr, g1, g2)
     if result == 0x02
         error("LibGEOS: Error in GEOSCoveredBy")
@@ -792,7 +792,7 @@ end
 # -----
 
 # GEOSGeometry ownership is retained by caller
-function prepareGeom(ptr::GEOSGeom, context::GEOSContext = _context)
+function prepareGeom(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSPrepare_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSPrepare")
@@ -800,7 +800,7 @@ function prepareGeom(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-function destroyPreparedGeom(ptr::Ptr{GEOSPreparedGeometry}, context::GEOSContext = _context)
+function destroyPreparedGeom(ptr::Ptr{GEOSPreparedGeometry}, context::GEOSContext = get_global_context())
     GC.@preserve context begin
         GEOSPreparedGeom_destroy_r(context.ptr, ptr)
     end
@@ -809,7 +809,7 @@ end
 function prepcontains(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedContains_r(context.ptr, g1, g2)
     if result == 0x02
@@ -821,7 +821,7 @@ end
 function prepcontainsproperly(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedContainsProperly_r(context.ptr, g1, g2)
     if result == 0x02
@@ -833,7 +833,7 @@ end
 function prepcoveredby(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedCoveredBy_r(context.ptr, g1, g2)
     if result == 0x02
@@ -845,7 +845,7 @@ end
 function prepcovers(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedCovers_r(context.ptr, g1, g2)
     if result == 0x02
@@ -857,7 +857,7 @@ end
 function prepcrosses(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedCrosses_r(context.ptr, g1, g2)
     if result == 0x02
@@ -869,7 +869,7 @@ end
 function prepdisjoint(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedDisjoint_r(context.ptr, g1, g2)
     if result == 0x02
@@ -881,7 +881,7 @@ end
 function prepintersects(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedIntersects_r(context.ptr, g1, g2)
     if result == 0x02
@@ -893,7 +893,7 @@ end
 function prepoverlaps(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedOverlaps_r(context.ptr, g1, g2)
     if result == 0x02
@@ -905,7 +905,7 @@ end
 function preptouches(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedTouches_r(context.ptr, g1, g2)
     if result == 0x02
@@ -917,7 +917,7 @@ end
 function prepwithin(
     g1::Ptr{GEOSPreparedGeometry},
     g2::GEOSGeom,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = GEOSPreparedWithin_r(context.ptr, g1, g2)
     if result == 0x02
@@ -939,7 +939,7 @@ end
 # -----
 # Unary predicate - return 2 on exception, 1 on true, 0 on false
 # -----
-function isEmpty(ptr::GEOSGeom, context::GEOSContext = _context)
+function isEmpty(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSisEmpty_r(context.ptr, ptr)
     if result == 0x02
         error("LibGEOS: Error in GEOSisEmpty")
@@ -947,7 +947,7 @@ function isEmpty(ptr::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function isSimple(ptr::GEOSGeom, context::GEOSContext = _context)
+function isSimple(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSisSimple_r(context.ptr, ptr)
     if result == 0x02
         error("LibGEOS: Error in GEOSisSimple")
@@ -955,7 +955,7 @@ function isSimple(ptr::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function isRing(ptr::GEOSGeom, context::GEOSContext = _context)
+function isRing(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSisRing_r(context.ptr, ptr)
     if result == 0x02
         error("LibGEOS: Error in GEOSisRing")
@@ -963,7 +963,7 @@ function isRing(ptr::GEOSGeom, context::GEOSContext = _context)
     result != 0x00
 end
 
-function hasZ(ptr::GEOSGeom, context::GEOSContext = _context)
+function hasZ(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSHasZ_r(context.ptr, ptr)
     if result == 0x02
         error("LibGEOS: Error in GEOSHasZ")
@@ -972,7 +972,7 @@ function hasZ(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Call only on LINESTRING (return 2 on exception, 1 on true, 0 on false)
-function isClosed(ptr::GEOSGeom, context::GEOSContext = _context)
+function isClosed(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSisClosed_r(context.ptr, ptr)
     if result == 0x02
         error("LibGEOS: Error in GEOSisClosed")
@@ -998,7 +998,7 @@ end
 #     GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE=1
 # };
 
-function isValid(ptr::GEOSGeom, context::GEOSContext = _context)
+function isValid(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSisValid_r(context.ptr, ptr)
     if result == 0x02
         error("LibGEOS: Error in GEOSisValid")
@@ -1028,7 +1028,7 @@ end
 # end
 
 # Return -1 on exception
-function geomTypeId(ptr::GEOSGeom, context::GEOSContext = _context)
+function geomTypeId(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeomTypeId_r(context.ptr, ptr)
     if result == -1
         error("LibGEOS: Error in GEOSGeomTypeId")
@@ -1037,7 +1037,7 @@ function geomTypeId(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Return 0 on exception
-function getSRID(ptr::GEOSGeom, context::GEOSContext = _context)
+function getSRID(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGetSRID_r(context.ptr, ptr)
     if result == 0
         error("LibGEOS: Error in GEOSGeomTypeId")
@@ -1045,7 +1045,7 @@ function getSRID(ptr::GEOSGeom, context::GEOSContext = _context)
     result
 end
 
-setSRID(ptr::GEOSGeom, context::GEOSContext = _context) = GEOSSetSRID_r(context.ptr, ptr)
+setSRID(ptr::GEOSGeom, context::GEOSContext = get_global_context()) = GEOSSetSRID_r(context.ptr, ptr)
 
 # May be called on all geometries in GEOS 3.x, returns -1 on error and 1
 # for non-multi geometries. Older GEOS versions only accept
@@ -1053,14 +1053,14 @@ setSRID(ptr::GEOSGeom, context::GEOSContext = _context) = GEOSSetSRID_r(context.
 # when fed simple geometries, so beware if you need compatibility with
 # old GEOS versions.
 """
-    numGeometries(geom, context=_context)
+    numGeometries(geom, context=get_global_context())
 
 Returns the number of sub-geometries immediately under a
 multi-geometry or collection or 1 for a simple geometry.
 For nested collections, remember to check if returned
 sub-geometries are **themselves** also collections.
 """
-function numGeometries(ptr::GEOSGeom, context::GEOSContext = _context)
+function numGeometries(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGetNumGeometries_r(context.ptr, ptr)
     if result == -1
         error("LibGEOS: Error in GEOSGeomTypeId")
@@ -1081,7 +1081,7 @@ Returns a copy of the specified sub-geometry of a collection.
 Numbering is one-based.
 For a simple geometry, returns a copy of the input.
 """
-function getGeometry(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
+function getGeometry(ptr::GEOSGeom, n::Integer, context::GEOSContext = get_global_context())
     n in 1:numGeometries(ptr, context) || error(
         "GEOSGetGeometryN: $(numGeometries(ptr, context)) sub-geometries in geom, therefore n should be in 1:$(numGeometries(ptr, context))",
     )
@@ -1098,12 +1098,12 @@ end
 Returns a vector of copy of the sub-geometries of a collection.
 For a simple geometry, returns the geometry in a vector of length one.
 """
-getGeometries(ptr::GEOSGeom, context::GEOSContext = _context) =
+getGeometries(ptr::GEOSGeom, context::GEOSContext = get_global_context()) =
     GEOSGeom[getGeometry(ptr, i, context) for i = 1:numGeometries(ptr, context)]
 
 # Converts Geometry to normal form (or canonical form).
 # Return -1 on exception, 0 otherwise.
-function normalize!(ptr::GEOSGeom, context::GEOSContext = _context)
+function normalize!(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSNormalize_r(context.ptr, ptr)
     if result == -1
         error("LibGEOS: Error in GEOSNormalize")
@@ -1112,7 +1112,7 @@ function normalize!(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Return -1 on exception
-function numInteriorRings(ptr::GEOSGeom, context::GEOSContext = _context)
+function numInteriorRings(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGetNumInteriorRings_r(context.ptr, ptr)
     if result == -1
         error("LibGEOS: Error in GEOSGetNumInteriorRings")
@@ -1121,7 +1121,7 @@ function numInteriorRings(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Call only on LINESTRING (returns -1 on exception)
-function numPoints(ptr::GEOSGeom, context::GEOSContext = _context)
+function numPoints(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeomGetNumPoints_r(context.ptr, ptr)
     if result == -1
         error("LibGEOS: Error in GEOSGeomGetNumPoints")
@@ -1130,7 +1130,7 @@ function numPoints(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Return -1 on exception, Geometry must be a Point.
-function getGeomX(ptr::GEOSGeom, context::GEOSContext = _context)
+function getGeomX(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     result = GEOSGeomGetX_r(context.ptr, ptr, out)
     if result == -1
@@ -1139,7 +1139,7 @@ function getGeomX(ptr::GEOSGeom, context::GEOSContext = _context)
     out[]
 end
 
-function getGeomY(ptr::GEOSGeom, context::GEOSContext = _context)
+function getGeomY(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     result = GEOSGeomGetY_r(context.ptr, ptr, out)
     if result == -1
@@ -1150,7 +1150,7 @@ end
 
 # Return NULL on exception, Geometry must be a Polygon.
 # Returned object is a pointer to internal storage: it must NOT be destroyed directly.
-function interiorRing(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
+function interiorRing(ptr::GEOSGeom, n::Integer, context::GEOSContext = get_global_context())
     if !(0 < n <= numInteriorRings(ptr, context))
         error(
             "LibGEOS: n=$n is out of bounds for Polygon with $(numInteriorRings(ptr, context)) interior ring(s)",
@@ -1163,7 +1163,7 @@ function interiorRing(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context
     cloneGeom(result, context)
 end
 
-function interiorRings(ptr::GEOSGeom, context::GEOSContext = _context)
+function interiorRings(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     n = numInteriorRings(ptr, context)
     if n == 0
         return GEOSGeom[]
@@ -1174,7 +1174,7 @@ end
 
 # Return NULL on exception, Geometry must be a Polygon.
 # Returned object is a pointer to internal storage: it must NOT be destroyed directly.
-function exteriorRing(ptr::GEOSGeom, context::GEOSContext = _context)
+function exteriorRing(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGetExteriorRing_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGetExteriorRing")
@@ -1183,7 +1183,7 @@ function exteriorRing(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Return -1 on exception
-function numCoordinates(ptr::GEOSGeom, context::GEOSContext = _context)
+function numCoordinates(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGetNumCoordinates_r(context.ptr, ptr)
     if result == -1
         error("LibGEOS: Error in GEOSGetNumCoordinates")
@@ -1192,7 +1192,7 @@ function numCoordinates(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
-function getCoordSeq(ptr::GEOSGeom, context::GEOSContext = _context)
+function getCoordSeq(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeom_getCoordSeq_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeom_getCoordSeq")
@@ -1202,19 +1202,19 @@ end
 # getGeomCoordinates(ptr::GEOSGeom) = getCoordinates(getCoordSeq(ptr))
 
 # Return 0 on exception (or empty geometry)
-getGeomDimensions(ptr::GEOSGeom, context::GEOSContext = _context) =
+getGeomDimensions(ptr::GEOSGeom, context::GEOSContext = get_global_context()) =
     GEOSGeom_getDimensions_r(context.ptr, ptr)
 
 # Return 2 or 3.
-getCoordinateDimension(ptr::GEOSGeom, context::GEOSContext = _context) =
+getCoordinateDimension(ptr::GEOSGeom, context::GEOSContext = get_global_context()) =
     GEOSGeom_getCoordinateDimension_r(context.ptr, ptr)
 
 """
-    getXMin(geom, context=_context)
+    getXMin(geom, context=get_global_context())
 
 Finds the minimum X value in the geometry
 """
-function getXMin(ptr::GEOSGeom, context::GEOSContext = _context)
+function getXMin(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     result = GEOSGeom_getXMin_r(context.ptr, ptr, out)
     if result == 0
@@ -1224,11 +1224,11 @@ function getXMin(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 """
-    getYMin(geom, context=_context)
+    getYMin(geom, context=get_global_context())
 
 Finds the minimum Y value in the geometry
 """
-function getYMin(ptr::GEOSGeom, context::GEOSContext = _context)
+function getYMin(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     result = GEOSGeom_getYMin_r(context.ptr, ptr, out)
     if result == 0
@@ -1238,11 +1238,11 @@ function getYMin(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 """
-    getXMax(geom, context=_context)
+    getXMax(geom, context=get_global_context())
 
 Finds the maximum X value in the geometry
 """
-function getXMax(ptr::GEOSGeom, context::GEOSContext = _context)
+function getXMax(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     result = GEOSGeom_getXMax_r(context.ptr, ptr, out)
     if result == 0
@@ -1252,11 +1252,11 @@ function getXMax(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 """
-    getYMax(geom, context=_context)
+    getYMax(geom, context=get_global_context())
 
 Finds the maximum Y value in the geometry
 """
-function getYMax(ptr::GEOSGeom, context::GEOSContext = _context)
+function getYMax(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     result = GEOSGeom_getYMax_r(context.ptr, ptr, out)
     if result == 0
@@ -1269,11 +1269,11 @@ end
 # TODO 02/2022: wait for libgeos release beyond 3.10.2 which will in include GEOSGeom_getExtent_r
 # Note: This is not in as of the current GEOS_jll 3.10.2...
 # """
-#     getExtent(geom, context=_context)
+#     getExtent(geom, context=get_global_context())
 
 # Retrieves the extent (xmin, ymin, xmax, ymax) of the geometry
 # """
-# function getExtent(ptr::GEOSGeom, context::GEOSContext = _context)
+# function getExtent(ptr::GEOSGeom, context::GEOSContext = get_global_context())
 #     out = Vector{Float64}(undef, 4)
 #     result = GEOSGeom_getExtent_r(
 #         context.ptr,
@@ -1290,7 +1290,7 @@ end
 # end
 
 # Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
-function getPoint(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
+function getPoint(ptr::GEOSGeom, n::Integer, context::GEOSContext = get_global_context())
     if !(0 < n <= numPoints(ptr, context))
         error(
             "LibGEOS: n=$n is out of bounds for LineString with numPoints=$(numPoints(ptr, context))",
@@ -1304,7 +1304,7 @@ function getPoint(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
 end
 
 # Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
-function startPoint(ptr::GEOSGeom, context::GEOSContext = _context)
+function startPoint(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeomGetStartPoint_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeomGetStartPoint")
@@ -1313,7 +1313,7 @@ function startPoint(ptr::GEOSGeom, context::GEOSContext = _context)
 end
 
 # Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
-function endPoint(ptr::GEOSGeom, context::GEOSContext = _context)
+function endPoint(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     result = GEOSGeomGetEndPoint_r(context.ptr, ptr)
     if result == C_NULL
         error("LibGEOS: Error in GEOSGeomGetEndPoint")
@@ -1324,7 +1324,7 @@ end
 # -----
 # Misc functions
 # -----
-function geomArea(ptr::GEOSGeom, context::GEOSContext = _context)
+function geomArea(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     # Return 0 on exception, 1 otherwise
     result = GEOSArea_r(context.ptr, ptr, out)
@@ -1334,7 +1334,7 @@ function geomArea(ptr::GEOSGeom, context::GEOSContext = _context)
     out[]
 end
 
-function geomLength(ptr::GEOSGeom, context::GEOSContext = _context)
+function geomLength(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     # Return 0 on exception, 1 otherwise
     result = GEOSLength_r(context.ptr, ptr, out)
@@ -1344,7 +1344,7 @@ function geomLength(ptr::GEOSGeom, context::GEOSContext = _context)
     out[]
 end
 
-function geomDistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function geomDistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     # Return 0 on exception, 1 otherwise
     result = GEOSDistance_r(context.ptr, g1, g2, out)
@@ -1354,7 +1354,7 @@ function geomDistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _contex
     out[]
 end
 
-function hausdorffdistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
+function hausdorffdistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context())
     out = Ref{Float64}()
     # Return 0 on exception, 1 otherwise
     result = GEOSHausdorffDistance_r(context.ptr, g1, g2, out)
@@ -1368,7 +1368,7 @@ function hausdorffdistance(
     g1::GEOSGeom,
     g2::GEOSGeom,
     densifyFrac::Real,
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     out = Ref{Float64}()
     # Return 0 on exception, 1 otherwise
@@ -1381,7 +1381,7 @@ end
 
 # Return 0 on exception, the closest points of the two geometries otherwise.
 # The first point comes from g1 geometry and the second point comes from g2.
-nearestPoints(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context) =
+nearestPoints(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = get_global_context()) =
     GEOSNearestPoints_r(context.ptr, g1, g2)
 
 # -----
@@ -1393,7 +1393,7 @@ nearestPoints(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context) =
 
 Return the size of the geometry's precision grid, 0 for FLOATING precision.
 """
-function getPrecision(geom::GEOSGeom, context::GEOSContext = _context)
+function getPrecision(geom::GEOSGeom, context::GEOSContext = get_global_context())
     GEOSGeom_getPrecision_r(context.ptr, geom)
 end
 
@@ -1421,7 +1421,7 @@ function setPrecision(
     geom::GEOSGeom,
     gridSize::Real,
     flags::Union{GEOSPrecisionRules,Integer},
-    context::GEOSContext = _context,
+    context::GEOSContext = get_global_context(),
 )
     result = geomFromGEOS(GEOSGeom_setPrecision_r(context.ptr, geom, gridSize, flags))
     if result == C_NULL
diff --git a/src/geos_operations.jl b/src/geos_operations.jl
index 625113c..1ff7e93 100644
--- a/src/geos_operations.jl
+++ b/src/geos_operations.jl
@@ -9,14 +9,14 @@ const GEOMTYPE = Dict{GEOSGeomTypes,Symbol}(
     GEOS_GEOMETRYCOLLECTION => :GeometryCollection,
 )
 
-writegeom(obj::Geometry, wktwriter::WKTWriter, context::GEOSContext = _context) =
+writegeom(obj::Geometry, wktwriter::WKTWriter, context::GEOSContext = get_global_context()) =
     _writegeom(obj.ptr, wktwriter, context)
-writegeom(obj::Geometry, wkbwriter::WKBWriter, context::GEOSContext = _context) =
+writegeom(obj::Geometry, wkbwriter::WKBWriter, context::GEOSContext = get_global_context()) =
     _writegeom(obj.ptr, wkbwriter, context)
-writegeom(obj::Geometry, context::GEOSContext = _context) = 
+writegeom(obj::Geometry, context::GEOSContext = get_global_context()) = 
     _writegeom(obj.ptr, context)
 
-function geomFromGEOS(ptr::GEOSGeom, context::GEOSContext = _context)
+function geomFromGEOS(ptr::GEOSGeom, context::GEOSContext = get_global_context())
     id = geomTypeId(ptr, context)
     if id == GEOS_POINT
         return Point(ptr, context)
@@ -38,33 +38,33 @@ function geomFromGEOS(ptr::GEOSGeom, context::GEOSContext = _context)
     end
 end
 
-readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = _context) =
+readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = get_global_context()) =
     geomFromGEOS(_readgeom(wktstring, wktreader, context), context)
-readgeom(wktstring::String, context::GEOSContext = _context) =
+readgeom(wktstring::String, context::GEOSContext = get_global_context()) =
     readgeom(wktstring, WKTReader(context), context)
 
-readgeom(wkbbuffer::Vector{Cuchar}, wkbreader::WKBReader, context::GEOSContext = _context) =
+readgeom(wkbbuffer::Vector{Cuchar}, wkbreader::WKBReader, context::GEOSContext = get_global_context()) =
     geomFromGEOS(_readgeom(wkbbuffer, wkbreader, context), context)
-readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = _context) =
+readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = get_global_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, context::GEOSContext = _context) =
+project(line::LineString, point::Point, context::GEOSContext = get_global_context()) =
     project(line.ptr, point.ptr, context)
-projectNormalized(line::LineString, point::Point, context::GEOSContext = _context) =
+projectNormalized(line::LineString, point::Point, context::GEOSContext = get_global_context()) =
     projectNormalized(line.ptr, point.ptr, context)
-interpolate(line::LineString, dist::Real, context::GEOSContext = _context) =
+interpolate(line::LineString, dist::Real, context::GEOSContext = get_global_context()) =
     Point(interpolate(line.ptr, dist, context), context)
-interpolateNormalized(line::LineString, dist::Real, context::GEOSContext = _context) =
+interpolateNormalized(line::LineString, dist::Real, context::GEOSContext = get_global_context()) =
     Point(interpolateNormalized(line.ptr, dist, context), context)
 
 # # -----
 # # Topology operations
 # # -----
 
-buffer(obj::Geometry, dist::Real, quadsegs::Integer = 8, context::GEOSContext = _context) =
+buffer(obj::Geometry, dist::Real, quadsegs::Integer = 8, context::GEOSContext = get_global_context()) =
     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 = _context,) =
+    context::GEOSContext = get_global_context(),) =
     geomFromGEOS(bufferWithStyle(obj.ptr, dist, quadsegs, endCapStyle, joinStyle, mitreLimit, context), context)
 
-envelope(obj::Geometry, context::GEOSContext = _context) =
+envelope(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(envelope(obj.ptr, context), context)
-envelope(obj::PreparedGeometry, context::GEOSContext = _context) =
+envelope(obj::PreparedGeometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(envelope(obj.ownedby.ptr, context), context)
-minimumRotatedRectangle(obj::Geometry, context::GEOSContext = _context) =
+minimumRotatedRectangle(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(minimumRotatedRectangle(obj.ptr, context), context)
-convexhull(obj::Geometry, context::GEOSContext = _context) =
+convexhull(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(convexhull(obj.ptr, context), context)
-boundary(obj::Geometry, context::GEOSContext = _context) =
+boundary(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(boundary(obj.ptr, context), context)
-unaryUnion(obj::Geometry, context::GEOSContext = _context) =
+unaryUnion(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(unaryUnion(obj.ptr, context), context)
-pointOnSurface(obj::Geometry, context::GEOSContext = _context) =
+pointOnSurface(obj::Geometry, context::GEOSContext = get_global_context()) =
     Point(pointOnSurface(obj.ptr, context), context)
-centroid(obj::Geometry, context::GEOSContext = _context) =
+centroid(obj::Geometry, context::GEOSContext = get_global_context()) =
     Point(centroid(obj.ptr, context), context)
-node(obj::Geometry, context::GEOSContext = _context) =
+node(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(node(obj.ptr, context), context)
 
-intersection(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+intersection(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(intersection(obj1.ptr, obj2.ptr, context), context)
-difference(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+difference(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(difference(obj1.ptr, obj2.ptr, context), context)
-symmetricDifference(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+symmetricDifference(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(symmetricDifference(obj1.ptr, obj2.ptr, context), context)
-union(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+union(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     geomFromGEOS(union(obj1.ptr, obj2.ptr, context), context)
 
 # # all arguments remain ownership of the caller (both Geometries and pointers)
@@ -123,52 +123,52 @@ union(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
 #     result
 # end
 
-simplify(obj::Geometry, tol::Real, context::GEOSContext = _context) =
+simplify(obj::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
     geomFromGEOS(simplify(obj.ptr, tol, context), context)
-topologyPreserveSimplify(obj::Geometry, tol::Real, context::GEOSContext = _context) =
+topologyPreserveSimplify(obj::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
     geomFromGEOS(topologyPreserveSimplify(obj.ptr, tol, context), context)
-uniquePoints(obj::Geometry, context::GEOSContext = _context) =
+uniquePoints(obj::Geometry, context::GEOSContext = get_global_context()) =
     MultiPoint(uniquePoints(obj.ptr, context), context)
-delaunayTriangulationEdges(obj::Geometry, tol::Real = 0.0, context::GEOSContext = _context) =
+delaunayTriangulationEdges(obj::Geometry, tol::Real = 0.0, context::GEOSContext = get_global_context()) =
     MultiLineString(delaunayTriangulation(obj.ptr, tol, true, context), context)
-delaunayTriangulation(obj::Geometry, tol::Real = 0.0, context::GEOSContext = _context) =
+delaunayTriangulation(obj::Geometry, tol::Real = 0.0, context::GEOSContext = get_global_context()) =
     GeometryCollection(delaunayTriangulation(obj.ptr, tol, false, context), context)
-constrainedDelaunayTriangulation(obj::Geometry, context::GEOSContext = _context) =
+constrainedDelaunayTriangulation(obj::Geometry, context::GEOSContext = get_global_context()) =
     GeometryCollection(constrainedDelaunayTriangulation(obj.ptr, context), context)
 
 
-sharedPaths(obj1::LineString, obj2::LineString, context::GEOSContext = _context) =
+sharedPaths(obj1::LineString, obj2::LineString, context::GEOSContext = get_global_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, context::GEOSContext = _context) =
+snap(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
     geomFromGEOS(snap(obj1.ptr, obj2.ptr, tol, context), context)
 
 # -----
 # Binary predicates
 # -----
 
-disjoint(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+disjoint(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     disjoint(obj1.ptr, obj2.ptr, context)
-touches(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+touches(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     touches(obj1.ptr, obj2.ptr, context)
-intersects(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+intersects(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     intersects(obj1.ptr, obj2.ptr, context)
-crosses(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+crosses(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     crosses(obj1.ptr, obj2.ptr, context)
-within(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+within(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     within(obj1.ptr, obj2.ptr, context)
-Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+Base.contains(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     Base.contains(obj1.ptr, obj2.ptr, context)
-overlaps(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+overlaps(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     overlaps(obj1.ptr, obj2.ptr, context)
-equals(obj1::Geometry, obj2::Geometry, context::GEOSContext = _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 = _context) =
+equalsexact(obj1::Geometry, obj2::Geometry, tol::Real, context::GEOSContext = get_global_context()) =
     equalsexact(obj1.ptr, obj2.ptr, tol, context)
-covers(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+covers(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     covers(obj1.ptr, obj2.ptr, context)
-coveredby(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+coveredby(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     coveredby(obj1.ptr, obj2.ptr, context)
 
 
@@ -176,27 +176,27 @@ coveredby(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
 # # Prepared Geometry Binary predicates
 # # -----
 
-prepareGeom(obj::Geometry, context::GEOSContext = _context) =
+prepareGeom(obj::Geometry, context::GEOSContext = get_global_context()) =
     PreparedGeometry(prepareGeom(obj.ptr, context), obj)
-Base.contains(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+Base.contains(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepcontains(obj1.ptr, obj2.ptr, context)
-containsproperly(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+containsproperly(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepcontainsproperly(obj1.ptr, obj2.ptr, context)
-coveredby(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+coveredby(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepcoveredby(obj1.ptr, obj2.ptr, context)
-covers(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context)= 
+covers(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context())= 
     prepcovers(obj1.ptr, obj2.ptr, context)
-crosses(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+crosses(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepcrosses(obj1.ptr, obj2.ptr, context)
-disjoint(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+disjoint(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepdisjoint(obj1.ptr, obj2.ptr, context)
-intersects(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+intersects(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepintersects(obj1.ptr, obj2.ptr, context)
-overlaps(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+overlaps(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepoverlaps(obj1.ptr, obj2.ptr, context)
-touches(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+touches(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     preptouches(obj1.ptr, obj2.ptr, context)
-within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context) =
+within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     prepwithin(obj1.ptr, obj2.ptr, context)
 
 # # -----
@@ -212,20 +212,20 @@ within(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = _context)
 # # -----
 # # Unary predicate - return 2 on exception, 1 on true, 0 on false
 # # -----
-isEmpty(obj::Geometry, context::GEOSContext = _context) =
+isEmpty(obj::Geometry, context::GEOSContext = get_global_context()) =
     isEmpty(obj.ptr, context)
-isEmpty(obj::PreparedGeometry, context::GEOSContext = _context) =
+isEmpty(obj::PreparedGeometry, context::GEOSContext = get_global_context()) =
     isEmpty(obj.ownedby.ptr, context)
-isSimple(obj::Geometry, context::GEOSContext = _context) =
+isSimple(obj::Geometry, context::GEOSContext = get_global_context()) =
     isSimple(obj.ptr, context)
-isRing(obj::Geometry, context::GEOSContext = _context) =
+isRing(obj::Geometry, context::GEOSContext = get_global_context()) =
     isRing(obj.ptr, context)
-isValid(obj::Geometry, context::GEOSContext = _context) =
+isValid(obj::Geometry, context::GEOSContext = get_global_context()) =
     isValid(obj.ptr, context)
-hasZ(obj::Geometry, context::GEOSContext = _context) =
+hasZ(obj::Geometry, context::GEOSContext = get_global_context()) =
     hasZ(obj.ptr, context)
 
-isClosed(obj::LineString, context::GEOSContext = _context) =
+isClosed(obj::LineString, context::GEOSContext = get_global_context()) =
     isClosed(obj.ptr, context) # Call only on LINESTRING
 
 # # -----
@@ -259,7 +259,7 @@ isClosed(obj::LineString, context::GEOSContext = _context) =
 # # -----
 
 # Gets the number of sub-geometries
-numGeometries(obj::Geometry, context::GEOSContext = _context) =
+numGeometries(obj::Geometry, context::GEOSContext = get_global_context()) =
     numGeometries(obj.ptr, context)
 
 # # Call only on GEOMETRYCOLLECTION or MULTI*
@@ -277,23 +277,23 @@ numGeometries(obj::Geometry, context::GEOSContext = _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 = _context) =
+getGeometry(obj::Geometry, n::Integer, context::GEOSContext = get_global_context()) =
     geomFromGEOS(getGeometry(obj.ptr, n, context), context)
-getGeometries(obj::Geometry, context::GEOSContext = _context) =
+getGeometries(obj::Geometry, context::GEOSContext = get_global_context()) =
     [geomFromGEOS(gptr, context) for gptr in getGeometries(obj.ptr, context)]
 
 # Converts Geometry to normal form (or canonical form).
-normalize!(obj::Geometry, context::GEOSContext = _context) =
+normalize!(obj::Geometry, context::GEOSContext = get_global_context()) =
     normalize!(obj.ptr, context)
 
 # LinearRings in Polygons
-numInteriorRings(obj::Polygon, context::GEOSContext = _context) =
+numInteriorRings(obj::Polygon, context::GEOSContext = get_global_context()) =
     numInteriorRings(obj.ptr, context)
-interiorRing(obj::Polygon, n::Integer, context::GEOSContext = _context) =
+interiorRing(obj::Polygon, n::Integer, context::GEOSContext = get_global_context()) =
     LinearRing(interiorRing(obj.ptr, n, context), context)
-interiorRings(obj::Polygon, context::GEOSContext = _context) =
+interiorRings(obj::Polygon, context::GEOSContext = get_global_context()) =
     map(LinearRing, interiorRings(obj.ptr, context))
-exteriorRing(obj::Polygon, context::GEOSContext = _context) =
+exteriorRing(obj::Polygon, context::GEOSContext = get_global_context()) =
     LinearRing(exteriorRing(obj.ptr, context), context)
 
 # # Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
@@ -321,33 +321,33 @@ exteriorRing(obj::Polygon, context::GEOSContext = _context) =
 #     result
 # end
 
-numPoints(obj::LineString, context::GEOSContext = _context) =
+numPoints(obj::LineString, context::GEOSContext = get_global_context()) =
     numPoints(obj.ptr, context) # Call only on LINESTRING
-startPoint(obj::LineString, context::GEOSContext = _context) =
+startPoint(obj::LineString, context::GEOSContext = get_global_context()) =
     Point(startPoint(obj.ptr, context), context) # Call only on LINESTRING
-endPoint(obj::LineString, context::GEOSContext = _context) =
+endPoint(obj::LineString, context::GEOSContext = get_global_context()) =
     Point(endPoint(obj.ptr, context), context) # Call only on LINESTRING
 
 # # -----
 # # Misc functions
 # # -----
 
-area(obj::Geometry, context::GEOSContext = _context) =
+area(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomArea(obj.ptr, context)
-geomLength(obj::Geometry, context::GEOSContext = _context) =
+geomLength(obj::Geometry, context::GEOSContext = get_global_context()) =
     geomLength(obj.ptr, context)
 
-distance(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) =
+distance(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) =
     geomDistance(obj1.ptr, obj2.ptr, context)
-hausdorffdistance(obj1::Geometry, obj2::Geometry, context::GEOSContext = _context) = 
+hausdorffdistance(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context()) = 
     hausdorffdistance(obj1.ptr, obj2.ptr, context)
 
-hausdorffdistance(obj1::Geometry, obj2::Geometry, densify::Real,context::GEOSContext = _context) =
+hausdorffdistance(obj1::Geometry, obj2::Geometry, densify::Real,context::GEOSContext = get_global_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, context::GEOSContext = _context)
+function nearestPoints(obj1::Geometry, obj2::Geometry, context::GEOSContext = get_global_context())
     points = nearestPoints(obj1.ptr, obj2.ptr, context)
     if points == C_NULL
         return Point[]
@@ -361,22 +361,22 @@ end
 # # Precision functions
 # # -----
 
-getPrecision(obj::Geometry, context::GEOSContext = _context) =
+getPrecision(obj::Geometry, context::GEOSContext = get_global_context()) =
     getPrecision(obj.ptr, context)
 setPrecision(
     obj::Geometry,
     grid::Real;
     flags = GEOS_PREC_VALID_OUTPUT,
-    context::GEOSContext = _context,) =
+    context::GEOSContext = get_global_context(),) =
     setPrecision(obj.ptr, grid::Real, flags, context)
 
 # ----
 #  Geometry information functions
 # ----
 
-getXMin(obj::Geometry, context::GEOSContext = _context) = getXMin(obj.ptr, context)
-getYMin(obj::Geometry, context::GEOSContext = _context) = getYMin(obj.ptr, context)
-getXMax(obj::Geometry, context::GEOSContext = _context) = getXMax(obj.ptr, context)
-getYMax(obj::Geometry, context::GEOSContext = _context) = getYMax(obj.ptr, context)
+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)
 
 # 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 e58bfbf..5f273f2 100644
--- a/src/geos_types.jl
+++ b/src/geos_types.jl
@@ -4,7 +4,7 @@ mutable struct Point <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a point from a pointer - only makese sense if it is a pointer to a point, otherwise error
-    function Point(ptr::GEOSGeom, context::GEOSContext = _context)
+    function Point(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         point = if id == GEOS_POINT
             point = new(cloneGeom(ptr, context), context)
@@ -16,11 +16,11 @@ mutable struct Point <: AbstractGeometry
         point
     end
     # create a point from a vector of floats
-    Point(coords::Vector{Float64}, context::GEOSContext = _context) =
+    Point(coords::Vector{Float64}, context::GEOSContext = get_global_context()) =
         Point(createPoint(coords, context), context)
-    Point(x::Real, y::Real, context::GEOSContext = _context) =
+    Point(x::Real, y::Real, context::GEOSContext = get_global_context()) =
         Point(createPoint(x, y, context), context)
-    Point(x::Real, y::Real, z::Real, context::GEOSContext = _context) =
+    Point(x::Real, y::Real, z::Real, context::GEOSContext = get_global_context()) =
         Point(createPoint(x, y, z, context), context)
 end
 
@@ -29,7 +29,7 @@ mutable struct MultiPoint <: AbstractGeometry
     context::GEOSContext
     # 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)
+    function MultiPoint(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         multipoint = if id == GEOS_MULTIPOINT
             new(cloneGeom(ptr, context), context)
@@ -47,7 +47,7 @@ mutable struct MultiPoint <: AbstractGeometry
         multipoint
     end
     # create a multipoint frome a vector of vector coordinates
-    MultiPoint(multipoint::Vector{Vector{Float64}}, context::GEOSContext = _context) =
+    MultiPoint(multipoint::Vector{Vector{Float64}}, context::GEOSContext = get_global_context()) =
         MultiPoint(
             createCollection(
                 GEOS_MULTIPOINT,
@@ -55,7 +55,7 @@ mutable struct MultiPoint <: AbstractGeometry
                 context),
             context)
     # create a multipoint from a list of points
-    MultiPoint(points::Vector{LibGEOS.Point}, context::GEOSContext = _context) =
+    MultiPoint(points::Vector{LibGEOS.Point}, context::GEOSContext = get_global_context()) =
         MultiPoint(
             createCollection(
                 GEOS_MULTIPOINT,
@@ -68,7 +68,7 @@ mutable struct LineString <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a linestring from a linestring pointer, otherwise error
-    function LineString(ptr::GEOSGeom, context::GEOSContext = _context)
+    function LineString(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         line = if id == GEOS_LINESTRING
             new(cloneGeom(ptr, context), context)
@@ -80,7 +80,7 @@ mutable struct LineString <: AbstractGeometry
         line
     end
     #create a linestring from a list of coordiantes
-    function LineString(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
+    function LineString(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context())
         line = new(createLineString(coords, context), context)
         finalizer(destroyGeom, line)
         line
@@ -91,7 +91,7 @@ mutable struct MultiLineString <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a multiline string from a multilinestring or a linestring pointer, else error
-    function MultiLineString(ptr::GEOSGeom, context::GEOSContext = _context)
+    function MultiLineString(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         multiline = if id == GEOS_MULTILINESTRING
             new(cloneGeom(ptr, context), context)
@@ -107,7 +107,7 @@ mutable struct MultiLineString <: AbstractGeometry
         multiline
     end
     # create a multilinestring from a list of linestring coordiantes
-    MultiLineString(multiline::Vector{Vector{Vector{Float64}}},context::GEOSContext = _context) =
+    MultiLineString(multiline::Vector{Vector{Vector{Float64}}},context::GEOSContext = get_global_context()) =
         MultiLineString(
             createCollection(
                 GEOS_MULTILINESTRING,
@@ -120,7 +120,7 @@ mutable struct LinearRing <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a linear ring from a linear ring pointer, otherwise error
-    function LinearRing(ptr::GEOSGeom, context::GEOSContext = _context)
+    function LinearRing(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         ring = if id == GEOS_LINEARRING
             new(cloneGeom(ptr, context), context)
@@ -133,7 +133,7 @@ mutable struct LinearRing <: AbstractGeometry
     end
     # create linear ring from a list of coordinates - 
     # first and last coordinates must be the same
-    function LinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
+    function LinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context())
         ring = new(createLinearRing(coords, context), context)
         finalizer(destroyGeom, ring)
         ring
@@ -145,7 +145,7 @@ mutable struct Polygon <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # 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, context::GEOSContext = _context)
+    function Polygon(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         polygon = if id == GEOS_POLYGON
             new(cloneGeom(ptr, context), context)
@@ -160,7 +160,7 @@ mutable struct Polygon <: AbstractGeometry
     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}}}, context::GEOSContext = _context)
+    function Polygon(coords::Vector{Vector{Vector{Float64}}}, context::GEOSContext = get_global_context())
         exterior = createLinearRing(coords[1], context)
         interiors = GEOSGeom[createLinearRing(lr, context) for lr in coords[2:end]]
         polygon = new(createPolygon(exterior, interiors, context), context)
@@ -168,10 +168,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 = _context) =
+    Polygon(ring::LinearRing, context::GEOSContext = get_global_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}, context::GEOSContext = _context) = 
+    Polygon(exterior::LinearRing, holes::Vector{LinearRing}, context::GEOSContext = get_global_context()) = 
         Polygon(
             createPolygon(exterior.ptr,
                           GEOSGeom[ring.ptr for ring in holes],
@@ -183,7 +183,7 @@ mutable struct MultiPolygon <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create multipolygon using a multipolygon or polygon pointer, else error
-    function MultiPolygon(ptr::GEOSGeom, context::GEOSContext = _context)
+    function MultiPolygon(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         multipolygon = if id == GEOS_MULTIPOLYGON
             new(cloneGeom(ptr, context), context)
@@ -203,7 +203,7 @@ mutable struct MultiPolygon <: AbstractGeometry
     end
 
     # create multipolygon from list of Polygon objects
-    MultiPolygon(polygons::Vector{Polygon}, context::GEOSContext = _context) =
+    MultiPolygon(polygons::Vector{Polygon}, context::GEOSContext = get_global_context()) =
         MultiPolygon(
             createCollection(
                 GEOS_MULTIPOLYGON,
@@ -212,7 +212,7 @@ mutable struct MultiPolygon <: AbstractGeometry
             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}}}}, context::GEOSContext = _context) =
+    MultiPolygon(multipolygon::Vector{Vector{Vector{Vector{Float64}}}}, context::GEOSContext = get_global_context()) =
         MultiPolygon(
             createCollection(
                 GEOS_MULTIPOLYGON,
@@ -230,7 +230,7 @@ mutable struct GeometryCollection <: AbstractGeometry
     ptr::GEOSGeom
     context::GEOSContext
     # create a geometric collection from a pointer to a geometric collection, else error
-    function GeometryCollection(ptr::GEOSGeom, context::GEOSContext = _context)
+    function GeometryCollection(ptr::GEOSGeom, context::GEOSContext = get_global_context())
         id = LibGEOS.geomTypeId(ptr, context)
         geometrycollection = if id == GEOS_GEOMETRYCOLLECTION
             new(cloneGeom(ptr, context), context)
@@ -242,7 +242,7 @@ mutable struct GeometryCollection <: AbstractGeometry
         geometrycollection
     end
     # create a geometric collection from a list of pointers to geometric objects
-    GeometryCollection(collection::Vector{GEOSGeom}, context::GEOSContext = _context) =
+    GeometryCollection(collection::Vector{GEOSGeom}, context::GEOSContext = get_global_context()) =
         GeometryCollection(
             createCollection(
                 GEOS_GEOMETRYCOLLECTION,
diff --git a/src/strtree.jl b/src/strtree.jl
index 0804cd6..49c545b 100644
--- a/src/strtree.jl
+++ b/src/strtree.jl
@@ -1,5 +1,5 @@
 """
-    STRtree(items; nodecapacity=10, context::GEOSContext=_context)
+    STRtree(items; nodecapacity=10, context::GEOSContext=get_global_context())
 
 Create a Sort Tile Recursive tree for fast intersection queries of objects for which an
 envelope is defined. The tree cannot be changed after its creation.
@@ -14,7 +14,7 @@ mutable struct STRtree{T}
     ptr::Ptr{GEOSSTRtree}  # void pointer to the tree
     items::T  # any geometry for which an envelope can be derived
     context::GEOSContext
-    function STRtree(items; nodecapacity = 10, context::GEOSContext = _context)
+    function STRtree(items; nodecapacity = 10, context::GEOSContext = get_global_context())
         tree = LibGEOS.GEOSSTRtree_create_r(context.ptr, nodecapacity)
         for item in items
             envptr = envelope(item).ptr
@@ -38,7 +38,7 @@ function destroySTRtree(obj::STRtree)
 end
 
 """
-    query(tree::STRtree, geometry; context::GEOSContext=_context)
+    query(tree::STRtree, geometry; context::GEOSContext=get_global_context())
 
 Returns the objects within `tree`, whose envelope intersects the envelope of `geometry`.
 


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index a50e918..f0b222c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,5 +21,6 @@ end
     include("test_regressions.jl")
     include("test_invalid_geometry.jl")
     include("test_strtree.jl")
+    include("test_misc.jl")
     Aqua.test_all(LibGEOS)
 end
diff --git a/test/test_geos_functions.jl b/test/test_geos_functions.jl
index 5d33bdb..c37e71b 100644
--- a/test/test_geos_functions.jl
+++ b/test/test_geos_functions.jl
@@ -8,7 +8,7 @@
 
     # round to 2 decimals
     writer = LibGEOS.WKTWriter(
-        LibGEOS._context,
+        LibGEOS.get_global_context(),
         trim = true,
         outputdim = 3,
         roundingprecision = 2,
@@ -17,7 +17,7 @@
 
     # round to 2 decimals and don't trim trailing zeros
     writer = LibGEOS.WKTWriter(
-        LibGEOS._context,
+        LibGEOS.get_global_context(),
         trim = false,
         outputdim = 3,
         roundingprecision = 2,
@@ -27,7 +27,7 @@
     # don't output the Z dimension
     p = readgeom("POINT(0.12345 2.000 0.1)")
     writer = LibGEOS.WKTWriter(
-        LibGEOS._context,
+        LibGEOS.get_global_context(),
         trim = false,
         outputdim = 2,
         roundingprecision = 2,
diff --git a/test/test_invalid_geometry.jl b/test/test_invalid_geometry.jl
index 4a2d4c8..a344cd2 100644
--- a/test/test_invalid_geometry.jl
+++ b/test/test_invalid_geometry.jl
@@ -9,6 +9,6 @@
     polygon = LibGEOS._readgeom("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),
         (15 15, 15 20, 20 20, 20 15, 15 15))")
     @test !LibGEOS.isValid(polygon)
-    @test LibGEOS.GEOSisValidReason_r(LibGEOS._context.ptr, polygon) ==
+    @test LibGEOS.GEOSisValidReason_r(LibGEOS.get_global_context().ptr, polygon) ==
           "Hole lies outside shell[15 15]"
 end
diff --git a/test/test_misc.jl b/test/test_misc.jl
new file mode 100644
index 0000000..2faf7a0
--- /dev/null
+++ b/test/test_misc.jl
@@ -0,0 +1,20 @@
+using Test
+using LibGEOS
+
+@testset "allow_global_context!" begin
+    Point(1,2,3)
+    LibGEOS.allow_global_context!(false)
+    @test_throws ErrorException Point(1,2,3)
+    Point(1,2,3, LibGEOS.GEOSContext())
+    LibGEOS.allow_global_context!(true) do
+        Point(1,2,3)
+    end
+    @test_throws ErrorException Point(1,2,3)
+    LibGEOS.allow_global_context!(true)
+    Point(1,2,3)
+    Point(1,2,3, LibGEOS.GEOSContext())
+    LibGEOS.allow_global_context!(false) do
+        @test_throws ErrorException Point(1,2,3)
+    end
+    Point(1,2,3)
+end