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@polygons[[1]]@Polygons, function(x) x@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@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: [email protected] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
