Polylabel.jl finds the pole of inaccessibility of a polygon, the most distant internal point from the polygon outline. This is useful for visual techniques like labelling polygons.
The main entry point is Polylabel.polylabel(polygon; atol, rtol) which processes any GeoInterface-compatible polygon (from GeometryBasics.jl, ArchGDAL.jl, LibGEOS.jl, Shapefile.jl, etc.) and returns a point as a 2-tuple of (x, y). It uses GeometryOps.jl to compute distances.
This algorithm was originally written (and taken from) mapbox/polylabel - you can find a lot more information there! To summarize, the algorithm is basically a quad-tree search across the polygon which finds the point which is most distant from any edge. There are alternative Julia implementations that are essentially the same algorithm in DelaunayTriangulation.jl
In the plot above, the pole of inaccessibility is shown in orange, while the input polygon (multipolygon in this case) is shown in blue.
First, get your polygon through whatever means:
using GeoInterface
p = GeoInterface.Polygon([[(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]])
# or load from a table from Shapefile, GeoJSON, GeoDataFrames, ArchGDAL, WellKnownGeometry, etc.Now, assuming p is your polygon or multipolygon (it can be from any GeoInterface package, like LibGEOS, ArchGDAL, GeometryBasics, Shapefile, GeoJSON, etc),
using Polylabel
polylabel(p) # (0.5, 0.5)will give you a result!
To shorten the time to compute, increase the keyword argument rtol (currently a 1% difference) or set atol to something at the scale you want.