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

[Karl Ove Hufthammer]

Karl Ove Hufthammer skreiv:

> 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

I now have an alternative, which takes into account the width and height 
of the label. It’s based on a brute-force search of candidate positions, 
and selects the position where the maximum distance between the label 
rectangle and the polygon edge is the greatest. Here’s an example of the 
function applied to a long and a tall string:

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

And here’s the source code:

-----
library(rgdal)
library(rgeos)

calc.labpt.strings=function(pol, label, gridpoints=2000) {
  # Fetch the label size
  wd=strwidth(label)
  ht=strheight(label)
  
  # Sample a regular grid of points inside the polygon
  co=coordinates(spsample(pol, n=gridpoints, type="regular"))
  
  # Create a rectangular polygon with a given size and at a given position
  makerect=function(x, y, wd, ht, rectID)
      Polygons(list(Polygon(cbind(c(x,x+wd,x+wd,x,x), c(y,y,y+ht,y+ht,y)))), rectID)
  
  # Create a candidate label rectangle for each grid point
  # (Note that ‘co’ may have fewer than ‘gridsize’ rows.)
  rects=SpatialPolygons(sapply(seq_len(nrow(co)), function(i) makerect(co[i,1], co[i,2], wd, ht, i)),
                        proj4string=CRS(proj4string(pol)))
  
  # Only keep the rectangles that are fully inside the polygon
  inside=apply(gContains(pol, rects, byid=TRUE), 1, any)
  if( all(!inside)) # Abort if no candidate label positions can fit the label
    stop("Could not fit label inside polygon (with the current number of gridpoints")
  rects=rects[inside,]
  
  # Convert the polygon to lines, and then measure the distance
  # from each label rectangle to the nearest line, keeping the
  # one with the largest distance.
  pol.l=as(pol, "SpatialLines")
  ind=which.max(apply(gDistance(pol.l, rects, byid=TRUE), 1, min))
  labelpos=apply(bbox(rects[ind,]), 1, mean)
  labelpos
}
-----


Here’s an example of how to use the function. Note that we have 
to draw the polygon *before* using the function, as the function
uses coordinate information from the graphical device.

The function uses 2000 candidate grid points by default, and is
still quite fast, but one can usually get by with fewer points
(even as few as 200 usually gets good results).

-----
# Fetch an example map
library(rworldmap)
kart=getMap(projection="equalArea")
xy.sp=kart["BRA",]

# Two examples
par(mfrow=c(1,2))

# A long (wide) string
plot(xy.sp, col="khaki")
label="A very long text string"
xy=calc.labpt.strings(xy.sp, label)
text(xy[1], xy[2], label)

# A tall string
plot(xy.sp, col="khaki")
label="N\na\nr\nr\no\nw"
xy=calc.labpt.strings(xy.sp, label)
text(xy[1], xy[2], label)
-----


Two final remarks:

The function seems to give excellent results on most *natural* 
maps, but for very regular polygons, my previous label 
placement function is better. This is because the optimality
criterion only looks at the maximum distance, and it doesn’t
handle ‘ties’ in a smart way either (it uses the first candidate 
point which satisfies the criterion).

The function will only work properly on projected polygons
(or polygons plotted with ‘asp=1’), and will warn if the
SpatialPolygons object is not projected. One good alternative 
for longlat data is to use the equirectangular projection, as
described in my previous post:

  http://article.gmane.org/gmane.comp.lang.r.geo/15170

I’d appreciate any comments and suggestions for improvements.

-- 
Karl Ove Hufthammer
E-mail: karl at huftis.org
Jabber: huftis at jabber.no