Generate centroids for vector polygons based on raster and land-use data using QGIS

Although QGIS has the ability to generate centroids for polygons, I am not sure how to generate centroids based on land-use development using raster data.

Here's an example of what I'm trying to do:

• I have a vector polygon layer, where each polygon represents a distinct zone for analysis.


• I have a color ortho raster layer that has the same geometric boundaries of my polygon layer.


• I am trying to develop a weighted centroid point layer based on the development density in the raster, yet still related somehow to the polygon it's originally in.



• Currently, I do this visually. It is time consuming to do 200+ zones and not as precise if a mathematical model was performing the analysis.


• I'm not sure if there is a plugin or a method for measuring the amount of roof space in a raster and using that and their density as a basis to develop centroids for the area

My Question:

How Can I generate centroids for vector polygons based on development density from a raster layer

Answer

When performing analyses of raster data polygon-by-polygon (with non-overlapping polygons), think first of zonal operations. By definition, these are statistical summaries for groups of cells identified by polygons, with an automatic iteration over the polygons.

OK, but exactly what summary? Look to the output for a clue: a centroid will be designated by an (x,y) coordinate pair. Therefore, you want some kind of weighted averages of coordinates.

Evidently, "development density" should serve as the weights. Coupling this with the definition of a weighted mean gives us a prototype algorithm:

1. 
   

   Create a grid of x-coordinates and a grid of y-coordinates. Call these [X] and [Y] respectively.

   
   


2. 
   

   Let [W] represent the weighting grid (development density). Compute [X.weighted] = [W] * [X] and [Y.weighted] = [W] * [Y]. (This is a local operation: it is performed cell by cell.)

   
   


3. 
   

   Obtain zonal sums of [X.weighted], [Y.weighted], and [W] itself. Instead of zonal sums, you may use zonal averages: the final results will not differ. (To do this, you may first have to rasterize the polygon layer.)

   
   


4. 
   
   

   Bearing in mind that the zonal sums will have one value for each polygon, join them on common polygon identifiers and (using a field calculation) divide the zonal sum of [X.weighted] by the zonal sum of [W] to obtain the weighted mean of [X]. Repeat, mutatis mutandis, for Y.

   
   

Now, QGIS uses GRASS GIS for raster operations. The GRASS zonal mean is performed by r.average. I suspect (but have not confirmed) that row() and col() in r.mapcalc will compute [X] and [Y] for you. Finally, r.mapcalc will also perform the local multiplications needed in Step (2).