[R-sig-Geo] Better label placement for polygons

[Karl Ove Hufthammer]

I have had some problems with bad label placements for SpatialPolygons,
so I wrote a small function for improving it. Hope it will be of use to
someone on this list. Here’s an example of its output:

  http://huftis.org/kritikk/polygon-labels.png

The orange points are the centroid values as returned by
‘coordinates()’ (note than one of them is outside its polygon), while
the red ones are the ones generated by the R function below.


The problem: The default label placements as returned by the
‘coordinates()’ function are the centroids. For convex polygons, they
work well, but for non-convex and non-symmetric polygons, the centroids
may be at a very different position from where one would manually place
the labels (and may even be *outside* the polygons).

A ‘good’ label position would be in the ‘bulk’ of the polygon, far away
from any edges. One way to find such a position is to use a negative
buffer to find interesting ‘inside’ areas of the polygon. This
automatically takes us away from any edges and closes any narrow
passages, which are the worst places for labels. My solution is thus to
use a buffer width so large that the *convex hull* of the resulting
interior polygon is wholly inside the original polygon. If the buffer
operation generates several polygons, I use the largest one (as this
would be the most natural place to put the label). The algorithm is
inspired by (but different from) the one is this article:
http://proceedings.esri.com/library/userconf/proc01/professional/papers/pap388/p388.html

Surprisingly, the resulting R code is quite fast (no doubt due to the
very efficient rgeos package!). For a simple polygon it takes just a
fraction of a second. And it can calculate label position for all the
polygons in the ‘rworldmap’ world polygons in a few seconds. The
resulting label placement also seems very good. Here’s the code:

-----
calc.labpt = function(pols) {
  # Prepopulate the label point matrix with the centroid values
  coords=coordinates(pols)
  
  # For each polygon in pols, calculate the appropriate label point
  for(i in seq_len(length(pols))) {
    
    # First fetch the polygon to process
    p=pols[i,]
    
    init=0                     # Initial amount to shrink
    estep=sqrt(gArea(p)/pi)/10 # Additional amount to shrink for each step
    
    # Try repeatedly shrinking the polygon until we’re left
    # with a polygon whose convex hulls fits inside 
    repeat {
      repeat {
        r = init + estep               # Amount to shrink
        p2 = gBuffer(p, width = -r)    # Shrink the polygon
        if( gArea(p2) <= 0 )           # If the shrunken polygon is empty ...
          estep = estep/2 else break   # ... try again with a smaller value
      }
      
      # If we’re left with more than one polygon, choose the largest one
      areas=sapply(p2 at polygons[[1]]@Polygons, function(x) x at area)
      if(length(areas) > 1) {
        # Note that we create a *new* SpatialPolygon containing the largest polygon.
        # I guess in theory we *could* have just replaced the @Polygons slot of p2,
        # but then gArea seems to crash R ... :(
        ind.max = which.max(areas)
        p2 = SpatialPolygons(list(Polygons(list(p2 at polygons[[1]]@Polygons[ind.max][[1]]),
                                         ID="middle")), proj4string=CRS(proj4string(p2)))
      }
      
      # Calculate the convex hull of the inner polygon.
      # If this is wholly contained in the original polygon,
      # break out of the loop and set the label point to
      # the centroid of the inner polygon.
      if( gContains(p, gConvexHull(p2)) ) break else init=init+estep
    }
    coords[i,] = coordinates(p2)
  }
  coords
}
-----

Here’s a simple example of how to use the function:

-----
library(rworldmap)
world = getMap(projection="equalArea")
norway = world["NOR",]

plot(norway)
points(coordinates(norway), col="orange", pch=19)
points(calc.labpt(norway), col="red", pch=19)
-----

Note 1: The function assumes that the polygons don’t contain any holes
(if they do, there may not exist a polygon calculated by buffering and
whose convex hull is wholly inside the original polygon). For most real
maps, where the holes are only small lakes, the function *may* very well
work (and usually will), but it’s not guaranteed to do so, and may enter
an endless loop.

Note 2: The label placement returned by the function is only (somewhat)
optimal *if* the label is shaped approximately like a square or
rectangle (e.g., a two-digit number). For very wide or tall labels,
there may exist much better label placements.

Note 3: The function assumes that the polygons are in a projected
coordinate system. When plotting unprojected polygons, the ‘plot()’
function automatically chooses an aspect ratio that make the result look
OK (as long as the latitude range is small). The label placement
function does not try to correct for this, and the result may look bad.
Example:

  norway = getMap(projection="none")["NOR",]
  plot(norway)        # ‘smart’ aspect ratio (looks OK)
  plot(norway, asp=1) # what ‘calc.labpt()’ sees
  
If you’re happy with the way ‘plot()’ displays your unprojected data,
just project it to the equivalent equirectangular projection before
calculating the label points:

  lat = mean(bbox(norway)[2,])                  # ‘Average’ latitude
  proj = CRS(paste0("+proj=eqc +lat_ts=", lat)) # Equirectangular
projection
  nor.proj = spTransform(norway, proj)          # Reproject polygon
  
‘plot(nor.proj)’ should *look* identical to ‘plot(nor)’ (but the
coordinates are very different). Now run ‘calc.labpt()’ on
‘nor.proj’ (and reproject back if you want to work in longlat
coordinates). (Note that for this example the result of ‘calc.labpt()’
on the original polygon is actually OK. In general, it will not be so.)

So, again, hope this will be useful for someone. It’s a quick hack, but
I think it works very well, and it has proven useful for me. I’d
appreciate any comments.

-- 
Karl Ove Hufthammer
http://huftis.org/
Jabber: karl at huftis.org