Calculating Weighted Centroids

By |

In this post, I will outline techniques for computing weighted-centroids in both QGIS and Google Earth Engine. For a polygon feature, the centroid is the geometric center. It can also be thought of as the average coordinate of all points within the polygon. There are some uses cases where you may want to compute a weighted-centroid where some parts of the polygon gets higher ‘weight’ than others. The main use-case is to calculate a population-weighted centroid. One can also use Night Lights data as a proxy for urbanized population and calculate a nightlights-weighted centroid. Some applications include:

  • Regional Planning: Locate the population-weighted centroid to know the most accessible location from the region.
  • Network Analysis: For generating demand points in location-allocation analysis, you need to convert demands from regions to points. It preferable to compute populated-weighted centroids for a more accurate analysis.

Do check out this twitter-thread by Raj Bhagat P for more discussion on weighted centroids.

Different Weighted Centroids for the State of Karnataka, India (2015)

Datasets

For this post, we will use the following datasets. If you want to follow along, download the data package and unzip the file on your computer.

Weighted Centroid Calculation in QGIS

QGIS Processing Toolbox comes with a built-in tool called Mean Coordinate that can be used to calculated weighted center from a point layer.

First, let’s calculate the Geometric Centroid. Load the karnataka layer and launch the Processing Toolbox → Centroids algorithm.

Save the result to a GeoPackage layer centroids.gpkg and name the layer geometric.

The resulting layer contains the centroid point of the state polygon.

Now, let’s compute the population-weighted centroid. Add the karnataka_ghsl_2015.tif to QGIS. This is a raster layer where pixel value is the number of people in the 1km x 1km grid.

We will convert the raster layer to a vector point layer. From the Processing Toolbox, run the Vector creation → Raster pixels to points algorithm.

Change the Field name to population and Run the tool.

A new point layer will be added to the Layers panel. There were no-data pixels in the source layer, so we have some points with an attribute value nan. We will filter them out. From the Processing Toolbox, launch the Vector selection → Extract by attribute algorithm.

Select ≠ as the Operator and enter nan as the Value. Click Run.

The resulting layer will have only those points with a valid value. Each point has an attribute called population. We will use that attribute to calculate the weighted centroid. Launch the Vector analysis → Mean Coordinate(s) algorithm.

Select population as the Weight field. This tool can also create multiple weighted-centroids – each calculated from a group of points. If you wanted to calculate centroids for each district in the state, and had an attribute with district name, you can choose that as the Unique ID field. For this exercise, leave it blank. Save the output to the centroids.gpkg with the layer name population.

The new layer population has the population-weighted centroid. You will notice that it is further towards the state capital which has a much higher share of the population.

You can repeat this process for the karnataka_viirs_2015 layer and calculate the Night Lights weighted centroid. This centroid is further towards the south, indicating more urbanized population towards the state capital.

Weighted Centroid Calculation in Google Earth Engine

The source datasets we used are global raster datasets. We can easily do this computation at regional or country level in Google Earth Engine. The The main function to be used is the ee.Image.pixelLonLat() which adds a new band to an image with the latitude and longitude of each pixel. This can be used to compute the weighted mean coordinates.

Below is the full script containing the code for pre-processing and weighted coordinate calculation. You can also access the Script in the Code Editor. Change the geometry variable to your region of interest to see the various centroids.

/*
Copyright (c) 2021 Ujaval Gandhi.
This work is licensed under the terms of the MIT license.  
For a copy, see https://opensource.org/licenses/MIT
*/
// Target Resolution
var outputScale = 5000;
// Target Projection
var outputProjection = ee.Projection('EPSG:4326').atScale(outputScale)
// Target Region
var admin1 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level1");
var filtered = admin1.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var geometry = filtered.geometry()
Map.centerObject(geometry)
// ****************** Pre-Process Population Data ********************
// *******************************************************************
var ghsl = ee.ImageCollection(
  "JRC/GHSL/P2016/POP_GPW_GLOBE_V1");
// Filter and get image for the year of interest
var filtered = ghsl
  .filter(ee.Filter.date('2015-01-01', '2016-01-01'))
var pop2015 = ee.Image(filtered.first())
// Aggregate to target resolution
var pop2015Scaled = pop2015
  .reduceResolution({
    reducer: ee.Reducer.sum().unweighted(),
    maxPixels: 1024
  })
  // Request the data at the scale and projection
  // of reduced resolution
  .reproject({
    crs: outputProjection
  });
var population = pop2015Scaled.rename('population')
// ****************** Pre-Process Night Lights Data ******************
// *******************************************************************
var nightlights = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG");
// Extract the projection before doing any computation
var projection = ee.Image(nightlights.first()).projection()
var filtered = nightlights.filter(ee.Filter.date('2015-01-01', '2016-01-01'))
var nightlights2015 = filtered.mean()
  .select('avg_rad').rename('nightlights')
// After we call mean() the projection is lost. Set it to source projection
var nightlights2015 = nightlights2015.setDefaultProjection(projection)
// Aggregate to target resolution
var nightlightsScaled = nightlights2015
  .reduceResolution({
    reducer: ee.Reducer.mean(),
    maxPixels: 1024
  })
  // Request the data at the scale and projection
  // of reduced resolution
  .reproject({
    crs: outputProjection
  });
var nightlights = nightlightsScaled.rename('nightlights')
// ****************** Clip and Visualize Images **********************
// *******************************************************************
var popClipped = population.clip(geometry)
var lightsClipped = nightlights.clip(geometry)
var popPalette = ['#fef0d9','#fdcc8a','#fc8d59','#e34a33','#b30000']
var lightsPalette = ['#253494','#2c7fb8', '#41b6c4', '#a1dab4', '#ffffcc']
Map.addLayer(popClipped, {min:500, max: 30000, palette: popPalette},
  'Population 2015', false)
Map.addLayer(lightsClipped, {min: 0.0, max: 60.0, palette: lightsPalette},
  'Night Lights 2015', false)
// ************************ Calculate Centroids **********************
// *******************************************************************
var image = popClipped.addBands(lightsClipped)
// Add lat and lon bands
var image = image.addBands(ee.Image.pixelLonLat())
// Sample Points
// The results will have properties extracted from each band of the image
var samples = image.sample({
  region: geometry,
  scale: outputScale,
  projection: outputProjection
})
// Calculate mean coordinates with different weight fields
var weightFields = ee.List(['nightlights', 'population'])
var centroids = weightFields.map(function(weightField){
    var xMean = ee.Number(samples.map(function(f) {
      var x = ee.Number(f.get('longitude'))
      var weight = ee.Number(f.get(weightField))
      return f.set('weightedx', x.multiply(weight))
    }).aggregate_sum('weightedx')).divide(
      samples.aggregate_sum(weightField))
  var yMean = ee.Number(samples.map(function(f) {
      var x = ee.Number(f.get('latitude'))
      var weight = ee.Number(f.get(weightField))
      return f.set('weightedy', x.multiply(weight))
    }).aggregate_sum('weightedy')).divide(
      samples.aggregate_sum(weightField))
  
  var centroid = ee.Geometry.Point([xMean, yMean])
  return ee.Feature(centroid, {'type': weightField})
})
var geometricCentroid = ee.Feature(
  geometry.centroid(10), {'type': 'geometric'})
var centroids = ee.FeatureCollection(
  centroids.add(geometricCentroid))
// Visualize Centroids
Map.addLayer(geometry, {color: 'gray'}, 'Region')
Map.addLayer(centroids, {color: 'red'}, 'Centroids')
Export.image.toDrive({
  image: popClipped,
  description: 'Population_Export',
  folder: 'earthengine',
  fileNamePrefix: 'karnataka_ghsl_2015',
  region: geometry,
  scale: outputScale,
  crs: outputProjection.crs().getInfo()
})
Export.image.toDrive({
  image: lightsClipped,
  description: 'Night_Lights_Export',
  folder: 'earthengine',
  fileNamePrefix: 'karnataka_viirs_2015',
  region: geometry,
  scale: outputScale,
  crs: outputProjection.crs().getInfo()
})
Export.table.toDrive({
  collection: centroids,
  description: 'Centroids_Export',
  folder: 'earthengine',
  fileNamePrefix: 'centroids',
  fileFormat: 'SHP'})

2 Comments

Leave a Comment

  1. Thank you so much for the explanation, Sir! I was following your steps to calculate population-weighted centroids for each of South Africa’s 52 districts. While it worked out for most of the district, unfortunately, there are two centroid points in some districts and none in others. Do you have any idea what the problem might be and how I could rectify it?

    Reply

Leave a Reply