Spatial statistics tools : clustering analysis on raster data

I've an apparently simple problem, but i can't find a clear methodology to use.

I'm tasked to delimit "urban areas" by vector convex polygons, using the Gridded Population of the World dataset from CIESIN

This dataset provides population density values over the entire world, as a raster file. The problem is, as you have already guessed, that the density values are very changing, and the definition of "urban" is quite relative.

I've tried to use a classical approach and calculated the slopes as if the density values were altitudes, but the slopes values were also very disparate and spatially complex, intricates.

I've looked into spatial clustering algorithms, LISA tools (Local Indocators of Spatial Association), with ArcGIS and GeoDa, but I'm quite lost among very specific tools. Some of the methods are working only on vector shapes, so a reclassification and a vectorization are needed (long computation).

Can you help me refining the set of methods and tools to use ? Thanks !

• What exactly is your definition of "urban"? BTW, I wouldn't expect all urban areas to be adequately described by convex polygons. Many have shapes controlled by geographic features--mountains, shorelines, and rivers--that are highly non-convex. – whuber Mar 20 '11 at 17:46
• The definition is not given, not a constant. I think the goal is to find substantically more dense areas in respect of their environment, which varies. So the idea about spatial statistics and LISA. You are right about the convexity, i should have wrote "not self-intersecting and not intersecting other polygons". – Laurent Jégou Mar 20 '11 at 19:05
• As a complement, i found a whole R module dedicated to spatial agglomeration tools : spdep. – Laurent Jégou Apr 7 '11 at 12:05

I did some work on this for my MSc http://ian01.geog.psu.edu/papers/mscthesis.pdf - basically I worked on gradient changes but the discussion may help you with this.

• Thanks, i'll read it soon, but that's promising :-) – Laurent Jégou Mar 20 '11 at 19:06
• I've tried the double Sobel 2nd derivative method you describe in the thesis, on a limited test zone, and it's very interesting ! Clear urban kernel are isolated, desptite very different density values. Many thanks :) By the way i used a free and opensource software to calculate the convolutions and raster math : Opticks. – Laurent Jégou Mar 22 '11 at 8:09
• A measure of artificial light would be a good proxy for urban conditions. A google search should find example studies. – b_dev May 3 '11 at 3:42
• @indiehacker - I looked at that but it's culturally dependant, France for example is much darker than the UK. – Ian Turton May 3 '11 at 13:36

From the point of view of population density, an "urban area" ought generally to satisfy just a few axiomatic criteria:

1. Its boundary should not include any points of (relatively) high density compared to the maximum density within its interior.

2. It should be simply connected (no "holes").

3. 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:

1. Optionally, smooth the density map.

2. Perform a "fill" algorithm on a density-related map (see below).

3. RegionGroup the result.

4. Remove holes from the RegionGrouped polygons.

5. Perform zonal sums of the population density over the filled polygons.

6. 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.

• Thanks for your detailed answer. I'll try to find the right software tools (or program them) to test the fill algorithm with the negative log, that seems a good lead. – Laurent Jégou Mar 22 '11 at 8:11