Thanks for sharing those details in the blog post!

Today I read that Mapbox also has an implementation of the "pole of inaccessibility" at https://github.com/mapbox/polylabel. And they also list other language implementations https://github.com/mapbox/polylabel#ports-to-other-languages.

A fast algorithm for finding the pole of inaccessibility of a polygon (in JavaScript and C++)

(Via: facebook/Rapid#534 (comment))

Since you referenced other implementation but not this one, I though I would share it.

@tordans The Mapbox implementation is the original one from Vladimir that I cite and link to in my blog post.

The computation worked well for my test area (Heidelberg city).

I like the approach with the ellipse a lot. In many cases it gives better results for labeling than the circle. Like in these examples:

However, in the following example I was surprised about the location of poi2. I would have expected it somewhere in the right-top of poi1 (blue)

However, in the following example I was surprised about the location of poi2. I would have expected it somewhere in the right-top of poi1 (blue)

If the polygon has several areas that are quite good it seems that it can take the algorithm a while to find the best one, it will be stuck on some local maximum somewhere and not consider that there is something better somewhere else. If you set the precision parameter smaller it should find the best position eventually but it will run much longer. You can try whether that makes a difference in this case.

This is available in PostGIS since 3.1.0: https://postgis.net/docs/ST_MaximumInscribedCircle.html

@pnorman had the suggestion in chat to scale by aspect ration of the bounding box (ST_YMax(g)-ST_YMin(g))/(ST_XMax(g)-ST_XMin(g)) instead of with a fixed scale.

Code here: https://gist.github.com/pnorman/58a58346915eb5c7369d8ff80bfac849

I've been playing with PIA via PostGIS in my latest style, and on the whole, it's better than other placement methods but does have some weaknesses.

My theory is that a better label placement algorithm is to take the inscribed ellipse with a maximum area. This can be done by rotating and scaling the polygon by different values, computing the maximum inscribed circle, undoing the rotation and scaling to compute the area, and then returning the point from the ellipse with the maximum area.

For practical reasons, it's not possible to try this for every possible scaling and rotation, so I've done it by scaling by 1.0, some hand-selected values, and the aspect ratio of the bounding box of the polygon. Assuming my theory above holds, this will return either the PIA or a point better than the PIA. I haven't tackled rotation, but my idea is to do a couple of rotations.

This is not perfect, but in my tests, does deal with the weaknesses of PIA while not causing any problems.

I'm not sure how we would handle the compute multiple points and select one based on properties of the computation in osm2pgsql.

I tried scaling with the the aspect ratio of the bounding box of the polygon and the result was sometimes better, sometimes worse. So for itself I don't think that's a win.

Calculating the PIA is expensive enough that I am sceptical about doing this many times for different aspect rations and/or rotations, so I don't think that's a good approach. I am certain there are better algorithms, but its unlikely we'll come up with something quickly. The main question for me is: Is the PIA we have here (with the scaling modification) a good enough solution for now that it makes sense baking it into osm2pgsql?

Calculating the PIA is expensive enough that I am sceptical about doing this many times for different aspect rations and/or rotations, so I don't think that's a good approach. I am certain there are better algorithms, but its unlikely we'll come up with something quickly. The main question for me is: Is the PIA we have here (with the scaling modification) a good enough solution for now that it makes sense baking it into osm2pgsql?

Yes. I'm not sure the scaling modification gains anything, but it's already coded.

I think this has been tested enough so I am removing the "work in progress" label. I just rebased this to the current master.

We can declare this feature as "experimental" in the docs, that gives us some room for changes in the future.

I have changed the Lua function so that it takes a Lua table with options. Currently only option is "stretch" (the second parameter from before). This allows us more flexibility when we want to change things later. The "precision" parameter is not exposed to the Lua code, because it is unclear how useful it is. Its really hard to explain what it does and users will probably just leave it at the default anyway. The now fixed (formally default) value is the same as with the PostGIS function.

I updated the docs in the example config.