Google Earth Engine makes it easy to compute statistics on gridded raster datasets. While calculating statistics on imagery datasets is easy, special care must be taken when working with population datasets. In this post, I will outline the correct technique for computing statistics for population rasters and aggregating pixels.

Calculating Total Population within a Region

The function reduceRegion() allows one to apply a statistical function (such as mean, median, count, sum etc.) on all pixels falling within the given region. You can use the ee.Reducer.sum() reducer to count the sum of all pixels values. Like many functions in Earth Engine, you must specify a scale value in the reduceRegion() function. The scale refers to the resolution at which the statistics will be computed. When computing sum on population datasets – always use the native resolution as the scale value. We can use the projection().nominalScale() function to find out the native resolution and use it for calculation.

var ghsl = ee.ImageCollection(
  "JRC/GHSL/P2016/POP_GPW_GLOBE_V1");

// Extract the projection before doing any computation
var projection = ee.Image(ghsl.first()).projection()
print('Native Resolution:', projection.nominalScale()) 
// 250 meters

// 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())

// Calculate Total Population for the City of Bengaluru
var bangalore = ee.FeatureCollection(
  "users/ujavalgandhi/public/bangalore_boundary");
var stats = pop2015.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: bangalore.geometry(),
  scale: projection.nominalScale(),
})
// Result is a dictionary
// Extract the value and print it
print(stats.get('population_count')) // 9557563.54557092

This is the correct value for the population in the region. But if you used a different scale value, you will get a very different and wrong result. Let’s understand why.

Scale and Image Pyramids

GIS users do not have to worry about scale or resolution when computing stats on rasters. Because all computation runs on the native scale. Earth Engine works slightly differently. To allow large computations, Earth Engine provides users with the flexibility of running computation at any scale – regardless of the resolution of the source image. When an image is ingested in Earth Engine, many lower resolution versions of the image are pre-computed. These are known as Image Pyramids. For each level of the pyramid, a 2×2 block of pixels is grouped into 1 larger pixel. When computing statistics at a specific scale, Earth Engine will use the appropriate pre-computed lower-resolution Pyramid Tiles to calculate the results.

The default pyramid tiles are created using mean value – lower resolution tiles are created with the average value of intersecting higher resolution pixels. This is called resampling and is ok for image datasets. But you cannot resample population datasets. For example, when you create a 1km resolution raster from a 500m resolution data – you want to sum the smaller pixel values to know the population within the larger pixel. This process is called aggregation.

Aggregating Pixels

To aggregate pixels in Earth Engine, you must use the reduceResolution() function. Learn more about this function in the Earth Engine User Guide. We will now see how to aggregate the population rasters into lower resolution rasters.

Why Aggregate Pixels?

There are many use cases where you would want to do this. Here are a few

• You want to compare different population datasets – each with a different pixel resolution. Aggregating them to a standard resolution will allow direct comparison between them. Read this article explaining the workflow.
• You want to overlay the population with other demographic datasets. Aggregating both datasets to the same resolution will allow you to do raster algebra operations.
• You want to compute stats for large regions. Aggregating pixels to a lower resolution will result in less number of pixels – reducing computation time.

We will now learn how to take the 250m resolution population raster from GHSL and aggregate it to 1km resolution.

First, we must compute the projection parameter at the required resolution. The projection consists of the CRS and a CRS transform. The transform consists of 6 parameters: [xScale, xShearing, xTranslation, yShearing, yScale, yTranslation]. We can use the helper function atScale() to calculate the new projection with the correct transform at the given scale.

Once we have the new projection, we can use the reduceResolution() function to aggregate the pixels using the given reducer function. When using the ee.Reducer.sum() reducer for aggregation, it is required to use an unweighted reducer ee.Reducer.sum().unweighted(). The reason for using an unweighted reducer here is due to the inner workings of Earth Engine as explained here and here.

// ****** Aggregation to 1km ****** //

// Get the projection at required scale
var projectionAt1k = projection.atScale(1000)

var pop2015At1k = pop2015
  .reduceResolution({
    reducer: ee.Reducer.sum().unweighted(),
    maxPixels: 1024
  })
  // Request the data at the scale and projection
  // of reduced resolution
  .reproject({
    crs: projectionAt1k
  });

var stats = pop2015At1k.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: geometry,
  scale: projectionAt1k.nominalScale(),
  maxPixels: 1e10,
  tileScale: 16})
print(stats.get('population_count')) // 9529826.09191158

Aggregation to Very Large Pixels

There is a limit of 1024 pixels in the reduceResolution function. That means a single step of aggregation can combine upto 1024 pixels in a single pixel. If the difference between source pixel size and target pixel size is large, you may get an error such as “Too many input pixels per output pixel“.

In such cases, you can do the aggregation in multiple steps. Here’s how to aggregate the 250m pixels image into a 10000m (10 km) pixel image.

// ****** Aggregation to 10km ****** //

// Get the projection at required scale
var projectionAt1k = projection.atScale(1000)
var projectionAt10k = projection.atScale(10000)

// Step1: 250m to 1000m
var scaledImage1 = pop2015
  .reduceResolution({
    reducer: ee.Reducer.sum().unweighted(),
    maxPixels: 1024
  })
  // Request the data at the scale and projection
  // of reduced resolution
  .reproject({
    crs: projectionAt1k
  });
// Step2: 1000m to 10000m
var pop2015At10k = scaledImage1
  .reduceResolution({
    reducer: ee.Reducer.sum().unweighted(),
    maxPixels: 1024
  })
  // Request the data at the scale and projection
  // of reduced resolution
  .reproject({
    crs: projectionAt10k
  });
  

Exporting Results

The aggregation operation can be computationally intensive and can time-out while calculating stats or displaying it on the map. The recommended practice is the export the result of the aggregation as an asset. Once the export finishes, you can import the asset into a script and continue to perform your computation. Earth Engine can also re-project the data to a specific CRS during the Export. Below is the code snippet showing how we can export the aggregated data to a UTM projection at 1km resolution.

// Clip and Export the results
// Export the aggregated data for a state
var admin1 = ee.FeatureCollection(
  "FAO/GAUL_SIMPLIFIED_500m/2015/level1");
var karnataka = admin1.filter(
  ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var geometry = karnataka.geometry()

// Source data is in World Mollweide projection
// Reproject it to a Local UTM projection for export
Export.image.toAsset({
  image: pop2015At1k.clip(geometry),
  description: 'Export_to_Asset',
  assetId: 'users/ujavalgandhi/public/karnataka_ghsl1k',
  region: geometry,
  scale: 1000,
  crs: 'EPSG:32643'
})
var exported = ee.Image(
  'users/ujavalgandhi/public/karnataka_ghsl1k')
Map.addLayer(exported,
  {min:500, max: 30000, palette: palette},
  'Exported Asset at 1km')

Here’s the full script with all the steps. If you have questions, do write them in the comments below.

If you are new to Earth Engine and want to master it, check out my course End-to-End Google Earth Engine.