Type instability of buffer function

[brunomarct]

Hello! First of all, thank you for creating and working on the JuliaGeo ecosystem. It is really awesome work!

The buffer function is not type stable even though LibGEOS.jl documentation indicates that it always returns a polygon and GEOS documentation specifies the same thing.

Edit: My initial understanding of the problem was not correct. See this comment below for the real question. Sorry for the messiness of the issue 😅

From my understanding:

• buffer(::Polygon, ::Real) returns a Polygon (edit: it can also return a MultiPolygon, see this comment below)
• buffer(::AbstractMultiGeometry, ::Real) returns a MultiPolygon (edit: it can also return a Polygon, see this comment below)

Here is two MWE:

For an AbstractGeometry:

mypoint = readgeom("POINT(1.0 1.0)")
buffer(mypoint, 1.0) |> typeof # Polygon

julia> @code_warntype buffer(mypoint, 1.0)
MethodInstance for buffer(::Point, ::Float64)
  from buffer(obj::Union{GeometryCollection, LineString, LinearRing, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}, dist::Real) @ LibGEOS ~/.julia/packages/LibGEOS/r6Xxl/src/geos_functions.jl:443
Arguments
  #self#::Core.Const(LibGEOS.buffer)
  obj::Point
  dist::Float64
Body::Union{GeometryCollection, LineString, LinearRing, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}
1 ─ %1 = LibGEOS.get_context(obj)::LibGEOS.GEOSContext
│   %2 = (#self#)(obj, dist, 8, %1)::Union{GeometryCollection, LineString, LinearRing, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}
└──      return %2

For an AbstractMultiGeometry:

mymultipoint = readgeom("MULTIPOINT(1.0 1.0, 5.0 5.0)")
buffer(mymultipoint, 1.0) |> typeof # MultiPolygon

julia> @code_warntype buffer(mymultipoint, 1.0)
MethodInstance for buffer(::MultiPoint, ::Float64)
  from buffer(obj::Union{GeometryCollection, LineString, LinearRing, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}, dist::Real) @ LibGEOS ~/.julia/packages/LibGEOS/r6Xxl/src/geos_functions.jl:443
Arguments
  #self#::Core.Const(LibGEOS.buffer)
  obj::MultiPoint
  dist::Float64
Body::Union{GeometryCollection, LineString, LinearRing, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}
1 ─ %1 = LibGEOS.get_context(obj)::LibGEOS.GEOSContext
│   %2 = (#self#)(obj, dist, 8, %1)::Union{GeometryCollection, LineString, LinearRing, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}
└──      return %2

I don't know if I missed any edge cases but I was wondering if there was anything that could be done to fix this?

[brunomarct]

Oh I just figured that buffer(::AbstractMultiGeometry, ::Real)does not always return MultiPolygon.

For example:

mymultipoint2 = readgeom("MULTIPOINT(1.0 1.0, 2.0 2.0)")
buffer(mymultipoint2, 1.0) |> typeof # Polygon

[brunomarct]

OK I just figured a case where buffer(::Polygon, ::Real)does not return a Polygon:

mypoly = readgeom("POLYGON ((0.5 2.5, 0.5 0.5, 4.5 0.5, 4.5 3.5, 3.5 3.5, 3.5 5.5, 4.5 5.5, 4.5 8.5, 0.5 8.5, 0.5 5.5, 2.5 5.5, 2.5 3.5, 0.5 3.5, 0.5 2.5))")
mypoly |> typeof # Polygon
buffer(mypoly, -0.5) |> typeof # MultiPolygon

So yeah, my previous understandings do not hold because of some edge cases.

Do you have any idea if there is a way to make buffer type stable for the special cases where its output is guaranteed to be a Polygon?

[asinghvi17]

We can always restrict output to Union{Polygon, MultiPolygon}, since operations like buffer(::Point, 0) still return empty polygons.

But it's not possible to know ahead of time whether the buffer of a geometry will return a Polygon or a MultiPolygon.

Is this causing any slowdowns for you? The memory representation of all these objects in Julia is the same (a pointer to the C/C++ GEOS object)...

[brunomarct]

Thanks for the rapid answer it answered my question.

I did not test whether this was causing any slowdowns. I'm still learning Julia, and when I was inspecting my code, I noticed the multiple possible return types of buffer. I thought this might indicate a problem or a potential speedup that could be gained. That led me to create this issue.

My initial use case was buffer(::Point, dist::Real) where dist > 0.0. I knew that the output would be a Polygon, but I wanted to make the compiler draw the same conclusion. With your answer and a bit more research, I've realized that there may not be a lot to gain, even if with the specialized code needed to achieve that.

Feel free to close the issue!