Calculating Area in Google Earth Engine

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/006/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.co.in/9c45ff677c46eae08952831de02bfb40

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.

23 Comments

Leave a Comment

    • To share your code, click the ‘Get Link’ button and paste the link. Also make any assets used in the code public so it can be accessed.

    • Hi, I have a similar question as yours. Have you solved your problem? Can you share with me how you fix the errors?

  1. This is excellent work, answering the question I had with the right mix of guidance, depth and complexity. I could feel that a deep knowledge of the subject meet with a great ability to teach.

  2. How can I calculate mean, mode, median, range, min, max of each band and export it into googleDrive in CSV format.

    Actually i am a beginner that’s why i don’t know anything about GEE. So Plz help me..

  3. Hi Ujaval,
    Excellent explanation and very explained.
    Can you also add a code snippet to calculate single class area but for all the districts (Area Calculation by Single Class by Admin Area.)
    Since I am absolutely new to both coding and GEE I am seeking this kind help from you.
    Kind Regards,
    Vikrant

  4. Thank you for the wonderful explanations on this post! I’ve learned a lot going through this blog and trying the steps on my own. I went to the End-to-end course you recommend at the conclusion of this blog and I’m working through it as well. I’m stuck at the 7th section though – the Direct Classification of Change. There are gcps of changed and unchanged pixels in the code. Are those generated by observation of the eye? Or there’s a step I’m to conduct between section 6 and 7 to obtain these?

    Thank you very much for any direction you can offer.

    Thanks to you and your team for making all these wonderful materials available. I’m learning so much. (And I’m on the wait-list for the live GEE course you offer. 🙂 )

    • Glad to know you found the materials useful. For the Direct Classification of Change, the points are generated by visual observation. These are the training points for the classifier which need to be generated where there is change and no-change. These can also come from field survey, but most people generate them by looking at the image pairs and finding pixels that changed or not. Hope that makes sense.

      • Oh yes, this is well understood. Thanks very much for taking time to respond.
        Best regards.

  5. There are no words to describe all the gratitude I feel for you spatialthoughts, THANK YOU Ujaval. You are the best.

  6. Hi ujaval, this is really awesome, I used your code (converted to python) last year and it worked well. Since yesterday, I am trying to compute the area on 2 classes and on a large image using the same script but I am getting the following error: ee.ee_exception.EEException: User memory limit exceeded. Here is my code

    # Class Area calculation

    def get_item(item):
    areaDict = ee.Dictionary(item)
    classNumber = ee.Number(areaDict.get(‘classification’)).format()
    # The result will be in square meters, this converts them into square kilometers
    area = ee.Number(areaDict.get(‘sum’)).divide(1e6).round()
    return ee.List([classNumber, area])

    def classAreaLists(class_areas):
    return class_areas.map(get_item)

    def export_class_area(classified):
    areaImage = ee.Image.pixelArea().addBands(classified)
    areas = areaImage.reduceRegion(**{
    “reducer”: ee.Reducer.sum().group(**{
    “groupField”: 1,
    “groupName”: “classification”,
    }),
    “geometry”: studysite.geometry(),
    “scale”: 10,
    “maxPixels”: 1e12
    })
    class_areas = ee.List(areas.get(“groups”))

    result = ee.Dictionary(classAreaLists(class_areas).flatten())
    print(result.getInfo())

    Do you have any idea on how I can solve this? so calling the getInfo() method on the dictionary crashes. I tried to increase the maxPixels but in vain.

    Thanks Hubert

    • Convert the result into a featurecollection and then Export it. Export can use more resources and your computation will finish

      Something like below
      fc = ee.FeatureCollection([ee.Feature(null, result)])

      See this for more infohttps://courses.spatialthoughts.com/end-to-end-gee.html#exporting-classification-results

  7. Hello dear Ujaval,

    I have been trying for days to find out why when I finish doing the supervised classification and I have applied the majorities filter and calculate the area of my class of interest, it gives different values after running the script.

    Also that generated me more doubts and I exported the result only of that class of interest and the value that ArcGIS calculates for example, is very different from the one obtained in GEE.

    What could be more or less the clue to solve this?

    I thank you in advance for any guidance you could give me.

    And thanks once again for everything.

    • Thank you very much Ujaval,

      I managed to solve it, what was happening is that I was comparing what it was not. I clarify that according to what I have noticed, the areas when compared with other software such as the one mentioned, these areas are usually slightly smaller.

      Have you noticed?

      • There are 2 reasons for the difference in the area.

        1. GEE uses a special method for area calculation using custom equal area projection for the region of the image – which is not what other software use. See https://groups.google.com/g/google-earth-engine-developers/c/Ccaorx-obVw/m/_ZQdP2wVAgAJ

        2. GEE can calculate area of fractional pixels. So if only 50% of a pixel is covered by your AOI, it will calculate it correctly. Whereas most other software (and data formats such as GeoTiff) have no way to represent fractional pixels and they count the number of whole pixels and multiply it by the number of pixels. This is where you’ll see slightly higher numbers of area.

Leave a Reply