0

I want to calculate the surfaces deforested every year for every protected area in Democratic Republic of the Congo. I am building from the script provided as answer by Rodrigo Principe in reply to this similar question. But this example used a private shapefile loaded in assets. When I try to reproduce this solution with a public dataset (in-land protected areas in DRC, so my question somehow completes the official GEE tutorial on quantigfying deforestation in DRC), it doesn't work : output "areas" remain desesperatelly empty. The code bellow is the same as proposed by Rodrigo Principe, except the lines marked with an"// EDIT" comment. Does anyone knows how to fix it?

// Load and filter the Hansen data
// EDIT : The source in following line has been updated to load v1_3
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015_v1_3')
              // END EDIT
              .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 = 100

//add country districts as a feature collection
// EDIT: loads public source+filters inland protected areas in DRC 
var distr = ee.FeatureCollection('WCMC/WDPA/current/polygons');
// Filter to keep only inland PA in 2 countries
var distr = distr.filter(ee.Filter.and(
    ee.Filter.eq("ISO3", "COD"),
    ee.Filter.eq("MARINE", "0")));
// END EDIT

//look at tree cover, find the area
var treeCover = gfc2014.select(['treecover2000']);
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())
               .divide(10000).select([0],["arealoss"]);

// total gain area
var gain = gfc2014.select(['gain'])
var areaGain = gain.gt(0).multiply(ee.Image.pixelArea())
               .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: 1e9 // EDIT : added to allow large areas
    })

    // 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")

print(areas)
  • best not to redefine "var distr" around line 17. delete 'var'. – intotecho Sep 20 '17 at 8:36
0

It does work, to prove it, I modified a little bit the code to show a chart. Just click on one feature and it shows 2 charts: loss and gain. If there are 2 overlapping features, it just plots for the first.

Code link: https://code.earthengine.google.com/7405deb2901a53729745cea43986e708

Code:

// Load and filter the Hansen data
// EDIT : The source in following line has been updated to load v1_3
var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015_v1_3')
              // END EDIT
              .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 = 100

//add country districts as a feature collection
// EDIT: loads public source+filters inland protected areas in DRC 
var distr = ee.FeatureCollection('WCMC/WDPA/current/polygons');
// Filter to keep only inland PA in 2 countries
var distr = distr.filter(ee.Filter.and(
    ee.Filter.eq("ISO3", "COD"),
    ee.Filter.eq("MARINE", "0")));
// END EDIT

//look at tree cover, find the area
var treeCover = gfc2014.select(['treecover2000']);
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())
               .divide(10000).select([0],["arealoss"]);

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

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

// Map cover area per feature
var districtSums = areaCover.reduceRegions({
  collection: distr,
  reducer: ee.Reducer.sum(),
  scale: scale,
});


var addVar = function(feature) {

  var nfeat = ee.Feature(ee.Feature(feature).geometry(), {})

  // function to iterate over the sequence of years
  var addVarYear = function(year, feat) {
    // cast var
    year = ee.Number(year).toInt()
    feat = ee.Feature(feat)
    //var nfeat = ee.Feature(feat.geometry(), {})

    // 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: 1e9 // EDIT : added to allow large areas
    })

    // 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 ee.Feature(feat.geometry(), {nemeloss:loss, namegain: gain})
    return feat.set(nameloss, loss, namegain, gain)
  }

  // iterate over the sequence
  var newfeat = ee.Feature(years.iterate(addVarYear, nfeat))

  // return feature with new properties
  return newfeat
}

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

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

var allY = ee.List.sequence(2001, 2014)
var allN = ee.List(['gain', 'loss'])
var allG = allY.map(function(item) {return ee.String('gain_').cat(ee.Number(item).int().format())})
var allL = allY.map(function(item) {return ee.String('loss_').cat(ee.Number(item).int().format())})

var onclick = function(coord) {
  var filt = areas.filterBounds(ee.Geometry.Point(coord.lon, coord.lat))
  var feat = ee.Feature(filt.first())
  var char_gain = ui.Chart.feature.byProperty({'features': feat,
                                               'xProperties': allG.getInfo()})
  print(char_gain.setOptions({'title':'GAIN'}))
  var char_loss = ui.Chart.feature.byProperty({'features': feat,
                                               'xProperties': allL.getInfo()})
  print(char_loss.setOptions({'title':'LOSS'}))
}
Map.onClick(onclick)
  • You are right ! The annual deforestation rates for every protected area polygon were there, it's just that I didn't understand how to access them. I'm not familiar with GeoJson featureCollections, which indeed handle data differently than shapefiles! And the point and click method you added to access the data is very nice. – fBedecarrats Sep 21 '17 at 22:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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