From the point of view of population density, an "urban area" ought generally to satisfy just a few axiomatic criteria:
Its boundary should not include any points of (relatively) high density compared to the maximum density within its interior.
It should be simply connected (no "holes").
Its average population density should exceed some prespecified threshold.
Axiom (1) is the most natural: if a boundary point were to have a high density, we would just move the boundary outward to include that point within the urban area. I would like to suggest that "relative" mean as a proportion of the maximum, such as one-tenth or one-hundredth or whatever. Axiom (2) avoids excluding parks and other low-density regions that naturally occur in cities. Axiom (3), which because it depends on a threshold is somewhat arbitrary, eliminates small compact villages.
Actually, there is at least one other element of arbitrariness: any gridded map of population density implicitly averages populations over local neighborhoods (equal to one cell in some cases and equal to a kernel radius for kernel density estimates). Let's accept this implicit neighborhood size (which can be changed by first running focal means or other kernel smooths over the original density map), this population threshold, and the sense of "relatively high" in axiom 1 as user-settable parameters that control the result.
These axioms lead naturally to a fairly simple algorithm: one must locate local maxima, look in their neighborhoods until a boundary is found to satisfy axiom (1), fill in any holes to satisfy axiom (2), and then screen all such candidate areas according to (3). This is done as follows:
Optionally, smooth the density map.
Perform a "fill" algorithm on a density-related map (see below).
RegionGroup the result.
Remove holes from the RegionGrouped polygons.
Perform zonal sums of the population density over the filled polygons.
Eliminate any polygons having sums (or mean densities) below the population (density) threshold.
What's left is your solution.
Let me say a little more about step (1), which is key. A fill algorithm identifies "sinks" and "fills" them up to a constant amount above their elevations. This is exactly what Axiom (1) asks us to do, provided (a) we can make "sink" play the role of "local maximum" and (b) make "constant amount above" play the role of "constant fraction of." The way to do this is by filling the negative logarithm of the density rather than the density itself. (Add a tiny constant first to the density--say, around 0.1 person per square kilometer--before taking the log, so that any cells containing zeros won't cause problems.) The "lakes" in the negative log density identify candidate urban areas. You still have three independent parameters to play with (entering at steps 0, 1, and 5); setting them will require some thought about what you really mean by "urban area" as well as some experimentation.
