Extracting building footprints from high-resolution imagery is a challenging task. Fortunately we now have access to ready-to-use building footprints dataset extracted using state-of-the-art ML techniques. Google’s Open Buildings project has mapped and extracted 1.8 billion buildings in Africa, South Asia, South-East Asia, Latin America and the Caribbean. Microsoft’s Global ML Building Footprints project has made over 1.24 billion building footprints available from most regions of the world.

Given the availability of these datasets, we can now analyze them to create derivative products. In this post, we will learn how to access these datasets and compute the aggregate count of buildings within a regular grid using Google Earth Engine. We will then export the grid as a shapefile and create a building density map in QGIS.

Here’s the overview of the workflow covered in this post

1. Select a region for mapping
2. Select the Buildings dataset (Google Open Buildings or Microsoft ML Buildings)
3. Create a regular grid of chosen size covering the region
4. Perform a Join to count the number of buildings in each grid
5. Export the resulting grid with building counts as a shapefile.
6. Apply symbology and create a map.

1. Select a Region for Mapping

You have many options on how to define your region of interest. In this post, we will select a Admin2 Region from the FAO GAUL dataset. Run the code snippet below in the GEE Code Editor.

var admin2 = ee.FeatureCollection(
  'FAO/GAUL_SIMPLIFIED_500m/2015/level2');
Map.addLayer(admin2);

Switch to the Inspector tab and click on a polygon. Expand the ‘▶ Feature …’ section and note the values for ADM0_NAME, ADM1_NAME and ADM2_NAME.

We will next filter the dataset to the chosen Admin2 region and extract it’s geometry.

var ADM0_NAME = 'India';
var ADM1_NAME = 'Karnataka';
var ADM2_NAME = 'Bangalore Urban';

var selected = admin2
  .filter(ee.Filter.eq('ADM0_NAME', ADM0_NAME))
  .filter(ee.Filter.eq('ADM1_NAME', ADM1_NAME))
  .filter(ee.Filter.eq('ADM2_NAME', ADM2_NAME))
var geometry = selected.geometry();
Map.addLayer(geometry, {color: 'gray'}, 'Selected Region');
Map.centerObject(geometry, 10);

If you do not want to use this method, you can simply drawing draw a polygon over your region using the Drawing Tools in the Code Editor instead of extracting it from an existing dataset. The script just expects a variable named geometry containing a polygon.

2. Select the Buildings Dataset

The Google Open Buildings dataset is available as a single FeatureCollection in the GEE Data Catalog. We can just load it and filter it for the chosen region. You can load the filtered layer to the map to check the coverage and quality. If your region is large, the layer may take a while to render. Zoom in closer to make the tiles load faster.

// Use the Google Open Buildings dataset
var googleBuildingsCol = ee.FeatureCollection(
  'GOOGLE/Research/open-buildings/v3/polygons');
var buildingsCol = googleBuildingsCol;
var buildings = buildingsCol.filter(ee.Filter.bounds(geometry));
Map.addLayer(buildings, {'color':'blue'}, 'Google Buildings');

The Microsoft Global Building Footprints dataset is available in the GEE Community Catalog. The dataset is divided into 1 FeatureCollection per country. So you will first need to query the asset folder to find the collection for your country. Once you find the collection, load it and then filter it for your chosen region.

var objects = ee.data.listAssets('projects/sat-io/open-datasets/MSBuildings')
print('Assets in MS Global Buildings Footprint Folder', objects['assets'])

// We are using the collection for India
var msBuildingsCol = ee.FeatureCollection(
  'projects/sat-io/open-datasets/MSBuildings/India');
var buildingsCol = msBuildingsCol;
var buildings = buildingsCol.filter(ee.Filter.bounds(geometry));
Map.addLayer(buildings, {'color':'red'}, 'MS Buildings');

3. Create a regular grid of chosen size covering the region

GEE API has a function coveringGrid() that can create a regular grid covering a polygon in a given CRS. Here we create a 1km x 1km grid in the EPSG:3857 CRS. This is a decent choice for most parts of the world, but you can adjust these values depending on your application.

// Creating a 1000m grid within the geometry
var gridSizeMeters = 1000;

// Generate a rectangular grid
var grid = geometry.coveringGrid('EPSG:3857', gridSizeMeters)

var gridOutline = ee.Image().byte().paint({
  featureCollection: grid,
  color: 1,
  width: 0.5
});
Map.addLayer(gridOutline, {palette: 'FF0000'}, 'Grids');

4. Perform a Join to count the number of buildings in each grid

We now need to count the number of buildings within each grid polygon. We can join the buildings layer with the grid layer based on intersecting features. We can then count the number of features within each grid.

var intersectFilter = ee.Filter.intersects({
    leftField: '.geo',
    rightField: '.geo',
    maxError: 10
});

var saveAllJoin = ee.Join.saveAll({
    matchesKey: 'buildings',
});

var joined = saveAllJoin
    .apply(grid, buildings, intersectFilter);
print(joined.first());

// Calculate total number of buildings within each feature.
var buildingCounts = joined.map(function(f) {
    var buildingsWithin = ee.List(f.get('buildings'));
    var totalBuildings = ee.FeatureCollection(buildingsWithin).size();
    return f.set('total_buildings', totalBuildings);
});

print(buildingCounts.first());

5. Export the resulting grid with building counts as a shapefile

Earth Engine let’s you run large computations interactively if they are able to finish within 5 minutes. For larger computations than this, you will have to Export the results. Filtering a collection containing billions of polygons and joining them with a grid is a fairly complex analysis and requires the resources of the batch compute to run. So we export the results. We are exporting a FeatureCollection and can choose any supported vector data format. Let’s export the grid with building counts as a shapefile. Run the code and start your Export task.

// We can export the grid with building count as a shapefile
// We use the 'selectors' option to only export
// the geometry and the total_buildings attribute
Export.table.toDrive({
  collection: buildingCounts.select(['total_buildings']),
  description: 'Grid_With_Building_Counts',
  folder: 'earthengine',
  fileNamePrefix: 'grid_with_building_counts',
  fileFormat: 'SHP',
})

6. Apply symbology and create a map

Once the export finishes, you will have a shapefile grid_with_building_counts.shp in your Google Drive. Download it to your computer. We can use any mapping software to create a map from this data. We will use QGIS. Open the shapefile in QGIS. You will see that the attribute table contains just a single column with building counts for each grid.

Use the Graduated renderer and apply a color ramp of your choice. You can define the class ranges manually based on the region. If you want a step-by-step tutorial for creating such maps – check out the Creating a Choropleth Map section of my Introduction to QGIS course.

Finally, you can create a Print Layout and create your map.

You can access the full script at https://code.earthengine.google.co.in/7e9bb75a9c0c49ce66b4b1e126a63123

If you end up using this script, do let me know in the comments.