Calculating Area in Google Earth Engine

By ujaval | June 19, 2020February 16, 2024

When working on Remote Sensing applications, many operations require calculating area. For example, one needs to calculate area covered by each class after supervised classification or find out how much area within a region is affected after a disaster. Calculating area for rasters and vectors is a straightforward operation in most software packages, but it is done in a slightly different way in Google Earth Engine – which can be confusing to beginners. In this post I will outline methods of calculating areas for both vectors as well as images. We will cover the following topics, starting from simple to complex.

Area Calculation for Features i.e. vector data
Area Calculation for Images (Single Class)
Area Calculation for Images by Class
Area Calculation for Images by Class by Region
Area Calculation for Images by Class by Region by Year

Preparing the Data

Land Cover Classification

For this post, we will use the MODIS 500m Landcover dataset from the Earth Engine Data Catalog. We will select the 2018 Land Cover Image from this collection. We pick the classification scheme in the band ‘LC_Type1’ which uses 17 different land use classes. The following code snippet shows how to select the image and visualize it.

var modisLandcover = ee.ImageCollection("MODIS/061/MCD12Q1")
var filtered = modisLandcover.filter(
  ee.Filter.date('2018-01-01', '2018-12-31'))
var landcover2018 = ee.Image(filtered.first())
var classified = landcover2018.select('LC_Type1')

var palette = ['05450a', '086a10', '54a708',
 '78d203', '009900', 'c6b044','dcd159', 
 'dade48', 'fbff13', 'b6ff05', '27ff87',
 'c24f44', 'a5a5a5', 'ff6d4c', '69fff8',
 'f9ffa4', '1c0dff']
Map.addLayer(classified,
 {min:1, max:17, palette: palette},
 'MODIS Landcover 2018')

Administrative Boundaries

To compute the area for regions, we will use the FAO GAUL: Global Administrative Unit Layers 2015, Second-Level Administrative Units dataset. This is a comprehensive global dataset that contains upto Level 2 (districts/counties/…) boundaries. For this post, we will filter it and use the district boundaries for the state of Kerala in India.

var gaul = ee.FeatureCollection(
  'FAO/GAUL_SIMPLIFIED_500m/2015/level2')
var kerala = gaul.filter(
  ee.Filter.eq('ADM1_NAME', 'Kerala'))
Map.addLayer(kerala, {color: 'purple'}, 'Admin2 Boundaries')

Now that we have polygons for the region of interest, we can clip the landcover image to the extent of this collection.

var keralaLandcover = classified.clip(kerala)
Map.addLayer(keralaLandcover,
 {min:1, max:17, palette: palette},
 'Kerala Land Cover 2018')

Area Calculation for Features

Calculating area for vector features is done using the built-in area() function. A pro-tip is that you can call .geometry() on a feature collection which gives the dissolved geometry of all features in the collection. We use this method to get a geometry representing the state of Kerala and compute its area.

var stateArea = kerala.geometry().area()
var stateAreaSqKm = ee.Number(stateArea).divide(1e6).round()
print(stateAreaSqKm)
// Result: 37999

Area Calculation for Images (Single Class)

Area calculation for images is done using the ee.Image.pixelArea() function. This function creates an image where each pixel’s value is the area of the pixel. If the image pixels contains values 0 or 1 – we can multiply this pixel area image with our image and calculate the total area using reduceRegion() function.

In this example, we extract the pixels for the class Urban Areas (class 13) and see how we can calculate the total built-up area in the state.

var urban = keralaLandcover.eq(13)
Map.addLayer(urban,
 {min:0, max:1, palette: ['grey', 'blue']},
 'Built-Up')

var areaImage = urban.multiply(ee.Image.pixelArea())
var area = areaImage.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: kerala.geometry(),
  scale: 500,
  maxPixels: 1e10
  })
var urbanAreaSqKm = ee.Number(
  area.get('LC_Type1')).divide(1e6).round()
print(urbanAreaSqKm)
// Result: 347

Area Calculation by Class

We learnt how to calculate area for a single class. But typically when you have a classified image, you want to compute area covered by each class. We must follow a similar process as before, but using a Grouped Reducer.

var areaImage = ee.Image.pixelArea().addBands(
      keralaLandcover)

var areas = areaImage.reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'class',
    }),
    geometry: kerala.geometry(),
    scale: 500,
    maxPixels: 1e10
    }); 

print(areas)

The result of reduceRegion() with a grouped reducer is a dictionary of dictionaries for each class. The top level dictionary has a single key named ‘groups‘. To extract the individual dictionaries and get properly formatted results, we must do a little post-processing. But before we dive into that, you must learn about another important function in Earth Engine called flatten().

flatten() is an important function in Earth Engine required for data processing. It takes a nested list and converts it to a flat list. Many Earth Engine constructors such a ee.Dictionary, ee.FeatureCollection etc. expect a flat list. So before creating such objects with nested objects, we must convert them to flat structures.

var nestedList = ee.List(
  [['a', 'b'], ['c', 'd'], ['e', 'f']])
print(nestedList) 
// Output: [["a","b"],["c","d"],["e","f"]]
print(nestedList.flatten())
// Output: ["a","b","c","d","e","f"]

Now we can take the results of the grouped reducer and map a function over it to extract the individual class areas and convert it to a single dictionary. Important to note that dictionary key must be of type ‘string’. Our keys are class numbers, so we use the format() method to convert the number to string

var classAreas = ee.List(areas.get('groups'))

var classAreaLists = classAreas.map(function(item) {
  var areaDict = ee.Dictionary(item)
  var classNumber = ee.Number(areaDict.get('class')).format()
  var area = ee.Number(
    areaDict.get('sum')).divide(1e6).round()
  return ee.List([classNumber, area])
})

var result = ee.Dictionary(classAreaLists.flatten())
print(result)

Area Calculation by Class by Region

We saw how we can calculate areas by class for the whole state. What if we wanted to know the breakup of these classes by each district? This requires one more level of processing. We can apply a similar computation as above, but wrap the computation in a function and by apply it using map() on the Feature Collection to obtain the values by each district geometry.

var calculateClassArea = function(feature) {
    var areas = ee.Image.pixelArea().addBands(classified)
    .reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'class',
    }),
    geometry: feature.geometry(),
    scale: 500,
    maxPixels: 1e10
    })

    var classAreas = ee.List(areas.get('groups'))
    var classAreaLists = classAreas.map(function(item) {
      var areaDict = ee.Dictionary(item)
      var classNumber = ee.Number(
        areaDict.get('class')).format()
      var area = ee.Number(
        areaDict.get('sum')).divide(1e6).round()
      return ee.List([classNumber, area])
    })

    var result = ee.Dictionary(classAreaLists.flatten())
    var district = feature.get('ADM2_NAME')
    return ee.Feature(
      feature.geometry(),
      result.set('district', district))
}

var districtAreas = kerala.map(calculateClassArea);

One thing to note is that each district may or may not have all of the 17 classes present. So each feature will have different number of attributes depending on which classes are present. We can explicitly set the expected fields in the output using the selectors argument for Export.table.toDrive() function. Because we need to use the list of output fields in the Export function we have to call getInfo() to get the list values on the client-side.

This result will be homogeneous table with all classes. Once done, we can export the results to a CSV file.

var classes = ee.List.sequence(1, 17)
var outputFields = ee.List(
    ['district']).cat(classes).getInfo()

Export.table.toDrive({
    collection: districtAreas,
    description: 'class_area_by_district',
    folder: 'earthengine',
    fileNamePrefix: 'class_area_by_district',
    fileFormat: 'CSV',
    selectors: outputFields
    })

Area Calculation by Class by Region by Year

If you have made it so far, you are ready to learn some more advanced techniques. What if you wanted to calculate the area for each district for multiple years? We can use one more map() function over a list of years, and run the computation for the whole MODIS Landcover collection.

We define a function that takes the year as the input and returns the district areas as the output. The code remains the same as before, but with one caveat. We add one more map() function and a ‘year’ property to each feature. This will help us identify the year of the calculation.

// Define a function that calculates class areas for 1 year
var calculateClassAreaByYear = function(year) {
  var startDate = ee.Date.fromYMD(year, 1, 1)
  var endDate = startDate.advance(1, 'year')
  var filtered = modisLandcover.filter(
    ee.Filter.date(startDate, endDate))
  
  var landcover = ee.Image(filtered.first())
  var classified = landcover.select('LC_Type1')
  
  // We write the map() function as before
  // But this time write it inline to make the code concise
  var districtAreas = kerala.map(function(feature) {
    
    var areas = ee.Image.pixelArea().addBands(classified)
      .reduceRegion({
        reducer: ee.Reducer.sum().group({
        groupField: 1,
        groupName: 'class',
      }),
      geometry: feature.geometry(),
      scale: 500,
      maxPixels: 1e10
      })
    var classAreas = ee.List(areas.get('groups'))
    var classAreaLists = classAreas.map(function(item) {
      var areaDict = ee.Dictionary(item)
      var classNumber = ee.Number(areaDict.get('class')).format()
      var area = ee.Number(areaDict.get('sum')).divide(1e6)//.round()
      return ee.List([classNumber, area])
    })
    var result = ee.Dictionary(classAreaLists.flatten())
    // The result dictionary has area for all the classes
    // We add the district name to the dictionary and create a feature
    
     // New! We need to now add the year in the property so we 
    // can identify which year did this calculation came from
    // Convert year to a string which looks better when we export.
    var district = feature.get('ADM2_NAME')
    return ee.Feature(feature.geometry(),
      result.set('district', district)
        .set('year', ee.Number(year).format('%d')))
  })
  return districtAreas
}

Now all that is left to do, is map() the above function on a list of years. Then we can export the results to a CSV.

// Now create a list of years and map() the above function
var years = ee.List.sequence(2001, 2018)
var results = years.map(calculateClassAreaByYear)

// Since there is a map()  within another map()
// We get a list of featurecollections, flatten it
// to create a featurecollection from the results.
var districtAreasByYear = ee.FeatureCollection(results).flatten()

// Set Output Fields, and make sure to include the 'year' field
var classes = ee.List.sequence(1, 17)

var outputFields = ee.List(['district', 'year']).cat(classes).getInfo()
// Export the results as a CSV file

Export.table.toDrive({
    collection: districtAreasByYear,
    description: 'class_area_by_district_by_year',
    folder: 'earthengine',
    fileNamePrefix: 'class_area_by_district_by_year',
    fileFormat: 'CSV',
    selectors: outputFields
    })

You can see the full script with comments at https://code.earthengine.google.com/6f0f7ed1aa57b5d52359b607d8910af6

Hope you found the post useful and got some inspiration to apply it to your own area calculation problem. If you are new to Earth Engine and want to master it, check out my course End-to-End Google Earth Engine.