1

I have a GEE-python code which first extracts the distinct pixel values and then aggregates the area of all such pixels within the adminitrative boundary. However, when running this function for big counties (e.g., Brazil and Argentina) the code shows computation limit error, probably due to unique_values = unique_values.getInfo().

When I try to do the same by looping the boundaries at adminitrative-level1 or level2, I still see the see the code taking longer. Is there an efficient way to achieve the same. I have been racking my brain for quite some time.

Code:

# Adding evaluated/processed 'area' bands to the Hansen dataset and calculating area
savefilename = 'NEW-PLANTATION-SP'
input = FL_att_to_new_plantations_species

variable = area_HansenLoss.where(input, area_HansenLoss);
Hansen_data_area = Hansen_data.addBands(variable)

# Calculate the sum of forest loss pixels for each feature in the geometry
area_geometrysum = variable.reduceRegions(
    collection=geometry,
    reducer=ee.Reducer.sum(),
    scale=Hansen_scale
)

species_dict = input.reduceRegion(
    reducer=ee.Reducer.frequencyHistogram(),
    geometry=geometry,
    scale=Hansen_scale,
    maxPixels=1e13
)

unique_values = ee.Dictionary(species_dict.get("species")).keys()

if unique_values.size().eq(0).getInfo():
    unique_values = ee.Array([0])
else:
    unique_values = ee.Array(unique_values.map(lambda val: ee.Number.parse(val)))

# Define the list of years to iterate over
startYear = 2001
endYear = int(str(datetime.datetime.now().year - 2))
years = ee.List.sequence(startYear - 2000, endYear - 2000)

# Function to add the forest loss data as properties to a feature
def addVar(feature, speciesclass):
    # Function to iterate over the sequence of years
    def addVarYear(year, feat):
        # Cast var
        year = ee.Number(year).toInt()
        feat = ee.Feature(feat)

        # Actual year to write as property
        actual_year = ee.Number(2000).add(year)

        # Filter the data by year
        filtered = Hansen_data_area.select("lossyear").eq(year).And(input.eq(speciesclass))

        # Apply the filter to the data
        filtered = Hansen_data_area.updateMask(filtered)

        # Reduce the forest loss data over the feature
        reduc = filtered.reduceRegion(
            geometry=feature.geometry(),
            reducer=ee.Reducer.sum(),
            scale=Hansen_scale,
            maxPixels=1e13
        )

        # Get the forest loss value
        loss = ee.Number(reduc.get("arealoss"))

        # Set the name of the property for the current year
        nameloss = ee.String('loss_').cat(actual_year.format())

        # Condition to handle cases where there is no forest loss data available
        cond = loss.gt(0)

        # Set the property for the current year
        return ee.Algorithms.If(
            cond,
            feat.set(nameloss, loss),
            feat
        )

    # Iterate over the sequence of years and add the properties to the feature
    newfeat = ee.Feature(years.iterate(addVarYear, feature))
    newfeat = newfeat.set("Species", speciesclass)
    
    # Return the feature with the new properties
    return newfeat

# Define an empty FeatureCollection to hold the output
areas = ee.FeatureCollection([])

unique_values = unique_values.getInfo()

# Loop over the speciesList and apply the addVar function with each speciesclass value
for speciesclass in unique_values:
    features = area_geometrysum.map(lambda feature: addVar(feature, speciesclass))
    areas = areas.merge(features)


# Print the resulting FeatureCollection
#print('FOREST LOSS TO ' + savefilename, outputFeatures);

# Generate the list of properties to export based on the start and end years
propertiesToExport = ['ADM0_CODE', 'ADM1_CODE', 'ADM2_CODE', 'Area_Boundary', 'Species']
for i in range(2001, endYear + 1):
    propertiesToExport.append('loss_' + str(i))

task = ee.batch.Export.table.toDrive(
            collection = areas,
            description = 'Forest_loss_to_'+savefilename +'_' + str(Admin_boundary) + '_' + '_' + str(datetime.datetime.now().isoformat()[0:19]),
            folder = 'Folder',
            fileFormat = 'CSV',
            selectors = propertiesToExport
            )
task.start()

The output looks like this (same code as above, but for small boundaries):enter image description here

1 Answer 1

2
+50

You can get rid of iteration by grouping when running reduceRegion(). In that way, after some data reformatting, you get area per loss year and species directly. You might still run in to trouble with large areas, but this solution is certainly more efficient than what you tried in your question.

regions = ee.FeatureCollection('FAO/GAUL_SIMPLIFIED_500m/2015/level2').filter(
    ee.Filter.eq('ADM0_NAME', 'Burundi')
)
lossYear = (
    ee.Image('UMD/hansen/global_forest_change_2021_v1_9').select('lossyear').add(2000)
)
species = ee.Image.random().multiply(3).round()  # Dummy species image


def to_region_features(region):
    def to_species_feature(species_group):
        species_group = ee.Feature(species_group)
        year_groups = ee.List(species_group.get('groups'))
        areas = ee.Array(
            year_groups.map(
                lambda properties: ee.Dictionary(properties).toArray().toList()
            )
        )
        labels = (
            areas.slice(1, 0, 1)
            .project([0])
            .toList()
            .map(lambda year: ee.String('loss_').cat(ee.Number(year).format('%d')))
        )
        values = areas.slice(1, 1).project([0]).toList()
        return (
            ee.Feature(
                None, region.toDictionary(['ADM0_CODE', 'ADM1_CODE', 'ADM2_CODE'])
            )
            .set('Species', species_group.get('Species'))
            .set(ee.Dictionary.fromLists(labels, values))
        )

    groups = ee.FeatureCollection(
        ee.List(
            ee.Image.pixelArea()
            .addBands(lossYear)
            .addBands(species)
            .reduceRegion(
                reducer=ee.Reducer.sum().group(1, 'lossYear').group(2, 'Species'),
                geometry=region.geometry(),
                scale=lossYear.projection().nominalScale(),
                maxPixels=1e13,
                tileScale=16,
            )
            .get('groups')
        ).map(
            # Include a size, to filter out cases where no data exist
            lambda group: ee.Feature(None, ee.Dictionary(group)).set(
                'size', 
                ee.List(ee.Dictionary(group).get('groups')).size()
            )
        )
    ).filter(ee.Filter.gt('size', 0))
    return ee.FeatureCollection(groups.map(to_species_feature))


areas = regions.map(to_region_features).flatten()
task = ee.batch.Export.table.toDrive(
    collection=areas,
    description='area-loss-by-region-and-species',
    selectors=['ADM0_CODE', 'ADM1_CODE', 'ADM2_CODE', 'Species'] + [
        'loss_{}'.format(lossYear) 
        for lossYear in range(2001, 2022)
    ],
).start()
3
  • I have two questions: (i) arealist is grouping the lossyear from Hansen, rather than grouping each year based on the plantation, right? (ii) If so, we still have to loop for unique values from the plantation dataset on the client side, which takes a similar processing time.
    – Ep1c1aN
    Apr 3, 2023 at 10:02
  • (i) Yes. (ii) Your script isn't executable, so it's hard to understand exactly what you want to do. You have an image with the species and you want the area by year and species? Apr 3, 2023 at 10:12
  • I updated my answer. Apr 3, 2023 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.