# Compute surfaces of intersection between circles in julia with LibGEOS

I need to compute surfaces of intersection between different circles. I am used to do it in R with the sf library and it works fine but I need to do it in Julia with the LibGEOS library and it is not as easy…

My first question would be, how to create a circle polygon? By creating the center and extend it with a buffer? For example:

My second question is, now that I have created some circles, can I put them in a common object, like a multipolygon?

My third and last question is, if I have an ensemble of circles, can I use the intersection function to compute all the intersections of all my circles?

I hope I am clear enough, I could really use some help regarding these questions!

Hi, welcome! Yeah this sure is possible. Not quite as well documented as sf though.

GEOS doesn’t have circles, only the geometry types listed here. Indeed taking a buffer of a point is a good way to create a circular polygon.

Why do you want to put the circles together in a common object? It should be doable to put them together in a MultiPolygon or GeometryCollection, but in the example below a simple Vector of Polygons is probably easier.

using LibGEOS, Plots

# 3 overlapping circles in a triangle
circles = [
buffer(Point(0.0, 0.0), 1.0),
buffer(Point(1.0, 0.0), 1.0),
buffer(Point(0.5, 1.0), 1.0),
]

isects = Polygon[]
for (i, g1) in enumerate(circles)
for (j, g2) in enumerate(circles)
# avoid counting double, and only use non empty intersections
if i > j && LibGEOS.intersects(g1, g2)
isect = LibGEOS.intersection(g1, g2)
push!(isects, isect)
end
end
end

plot(isects)

I guess for the intersection function I wanted an object like a multipolygon because I need to compute all the intersections between the circles, not only 2 by 2 intersections. With your example, I would also need to have a polygon for the intersection between the three circles. Ideally I would need as many polygons as overlapped surfaces, just like this figure below:

But since the intersection function only takes 2 polygons as arguments it might not be possible?

Hmm, I’m actually not sure how to do this. So in your figure there are no overlapping areas anymore right? It looks like the result of a union in QGIS. However union in GEOS also directly dissolves the result. Though probably with a combination of the other functions you can effectively do the same?

# need to go through GeoInterface, see issue #87
mp = MultiPolygon(GeoInterface.coordinates(circles))
plot(LibGEOS.unaryUnion(mp))

You should be able to write each colored surface as intersections between circles or their complements.
For example, the orange region is A \cap B \cap \neg C, where A is the top circle, B is the right circle and C is the left circle.

1 Like