Sunday 25 August 2019

Iterating over years for features in feature collection using Google Earth Engine?


I am currently working with the Hansen data (Global Forest Change) in Earth Engine. I have also imported a fusion table representing districts in a specific country. My goal is as follows:



I would like to create a table which aggregates deforestation data for each district in the feature table. For example, a given district ('feature') should eventually have the following associated properties: sum of gained forest area for each year, sum of lost forest area for each year, and base year tree cover. Once all of those are gathered, I plan to export the data.


I have figured out how to do this for one year, but the problem is I can't figure out how to iterate over years without doing everything manually. I am also not entirely sure whether my filtering actually gives me the year that I expect, as I copied the filtering from an example in the tutorials.


Here is some of my current code, which aggregates the tree cover for each district, and adds a 'loss' property to that feature collection:


//Load and filter the Hansen data
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015').select(['treecover2000','loss','gain','lossyear']);

//add country districts as a feature collection
var distr = ee.FeatureCollection('ft:1U7sXFHXtxQ--g7XMeXlvPhNXPBcDtPg8Yzr2pvsg', 'geometry');

//look at tree cover, find the area

var treeCover = gfc2014.select(['treecover2000']);
var areaCover = treeCover.multiply(ee.Image.pixelArea());

var lossIn2001 = gfc2014.select(['lossyear']).eq(1);
var areaLoss = lossIn2001.multiply(ee.Image.pixelArea());

//need to check that .eq(1) is actually returning data for
//2001
var gainIn2001 = gfc2014.select(['gain']).eq(1);
var areaGain = gainIn2001.multiply(ee.Image.pixelArea());


var districtSums = areaCover.reduceRegions({
collection: distr,
reducer: ee.Reducer.sum(),
scale: 100,
});

//This function computes pixels lost for each district
var addLoss = function(feature) {
var loss = areaLoss.reduceRegion({

reducer: ee.Reducer.sum(),
geometry: feature.geometry(),
scale: 100
});
return feature.set({loss2001: loss.get('lossyear')});
};

// Map the area getting function over the FeatureCollection.
var areaLosses = districtSums.map(addLoss);


I have looked to the following pages for help, but unfortunately neither appears to answer my question:



This is my first time working with geographic data so the data structures are not intuitive to me yet.



Answer



You are filtering year in the correct way. This is how I'd do it:


//Load and filter the Hansen data
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015')
.select(['treecover2000','loss','gain','lossyear']);

// list for filter iteration

var years = ee.List.sequence(1, 14)

// turn your scale into a var in case you want to change it
var scale = gfc2014.projection().nominalScale()

//add country districts as a feature collection
var distr = ee.FeatureCollection('ft:1U7sXFHXtxQ--g7XMeXlvPhNXPBcDtPg8Yzr2pvsg', 'geometry');

//look at tree cover, find the area
var treeCover = gfc2014.select(['treecover2000']);


// most recent version of Hansen's data has the treecover2000 layer
// ranging from 0-100. It needs to be divided by 100 if ones wants
// to calculate the areas in ha and not hundreds of ha. If not, the
// layers areaLoss/areaGain are not comparable to the areaCover. Thus
treeCover = treeCover.divide(100); // Thanks to Bruno

var areaCover = treeCover.multiply(ee.Image.pixelArea())
.divide(10000).select([0],["areacover"])


// total loss area
var loss = gfc2014.select(['loss']);
var areaLoss = loss.gt(0).multiply(ee.Image.pixelArea()).multiply(treeCover)
.divide(10000).select([0],["arealoss"]);

// total gain area
var gain = gfc2014.select(['gain'])
var areaGain = gain.gt(0).multiply(ee.Image.pixelArea()).multiply(treeCover)
.divide(10000).select([0],["areagain"]);


// final image
var total = gfc2014.addBands(areaCover)
.addBands(areaLoss)
.addBands(areaGain)

Map.addLayer(total,{},"total")

// Map cover area per feature
var districtSums = areaCover.reduceRegions({
collection: distr,

reducer: ee.Reducer.sum(),
scale: scale,
});


var addVar = function(feature) {

// function to iterate over the sequence of years
var addVarYear = function(year, feat) {
// cast var

year = ee.Number(year).toInt()
feat = ee.Feature(feat)

// actual year to write as property
var actual_year = ee.Number(2000).add(year)

// filter year:
// 1st: get mask
var filtered = total.select("lossyear").eq(year)
// 2nd: apply mask

filtered = total.updateMask(filtered)

// reduce variables over the feature
var reduc = filtered.reduceRegion({
geometry: feature.geometry(),
reducer: ee.Reducer.sum(),
scale: scale,
maxPixels: 1e13
})


// get results
var loss = ee.Number(reduc.get("arealoss"))
var gain = ee.Number(reduc.get("areagain"))

// set names
var nameloss = ee.String("loss_").cat(actual_year)
var namegain = ee.String("gain_").cat(actual_year)

// alternative 1: set property only if change greater than 0
var cond = loss.gt(0).or(gain.gt(0))

return ee.Algorithms.If(cond,
feat.set(nameloss, loss, namegain, gain),
feat)

// alternative 2: always set property
// set properties to the feature
// return feat.set(nameloss, loss, namegain, gain)
}

// iterate over the sequence

var newfeat = ee.Feature(years.iterate(addVarYear, feature))

// return feature with new properties
return newfeat
}

// Map over the FeatureCollection
var areas = districtSums.map(addVar);

Map.addLayer(areas, {}, "areas")


In that script you get 3 fields: loss_{year}, gain_{year}, sum But if you want better 4 fields: loss, gain, year, sum; change for:


return ee.Algorithms.If(cond, 
feat.set("loss", loss, "gain", gain, "year", actual_year),
feat)

You could also compute percentage and set it to the features.


Edit: Thank to @Bruno_Conte_Leite, who made me reconsider my answer, I have made some updates, the one suggested by Bruno and others.





  1. Scale: I suggest to keep the original scale of Hansen data.




  2. treeCover: most recent version of Hansen's data has the treecover2000 layer ranging from 0-100. It needs to be divided by 100 if ones wants to calculate the areas in ha and not hundreds of ha. (Bruno)




  3. areaLoss and areaGain: Added .multiply(treeCover) otherwise the area would be of the whole pixel and not of the indicated percentage




  4. maxPixels: I added maxPixels: 1e13 in the reduction





No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...