1

I have a working script that calculates the number of times a 30 m pixel intersects with a polygon feature in a ee.FeatureCollection. The input data is a 30 m raster for the entire mainland SE Asian region and a polygon layer. To generate the number of polygon counts per pixel, I create an Image Collection where each image corresponds to a raster with values of '1' for each pixel within a single polygon feature, and after ensuring that all the images share the same image footprint, I sum all the values in the ee.ImageCollection.

The script works (link here; I shared the assets too), and I have managed to get the result for Cambodia, which was an image of size 10 Mb.

//  Insert MSEA country boundaries
var msea = ee.FeatureCollection('users/jjohanness1992/GADM/mainland_SEA');
var country = msea.filter(ee.Filter.eq('ISO','LAO'));
print('country metadata',country);

//  Load image
var image = ee.Image('users/jjohanness1992/Fire_and_LCC/clip_lcc_2011_2012_templateLCC').clip(country);
// this is a 30 m raster for the entire mainland SE Asia

//  Create reclassification codes
//  Here I reclassify the land cover map to a raster of 0's
//  This is because I want to preserve this raster grid when performing the polygon count
var oldgroup = ee.List([18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306
]);
var newgroup = ee.List([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
]);

//  Reclassify lcc image to an image of 0 values
var image2 = image.remap(oldgroup,newgroup).select('remapped').rename('b1').toInt();


//  Load polygons
var polygons = ee.FeatureCollection('users/jjohanness1992/Fire_and_LCC/msea_fire_pixels_2012').filterBounds(country); // contains 3 fire polygons

//  Define function to iterate per fire polygon feature
var firepolygoncountperpixel = function(feature) {
  var f = feature;
  var i = image2.clip(f).toInt();
  var r = i.remap({from: ee.List([0]), to: ee.List([1])}).select('remapped').rename('b1').toInt();
  var c = ee.ImageCollection.fromImages([image2, r]);
  var m = c.mosaic().toInt();
  return m;
};

//  Apply function over all features
var mosaic_coll = polygons.map(firepolygoncountperpixel);
//print('mosaic_coll',mosaic_coll);
//  Result is an ImageColl comprised of 'x' images with 1 band each

//  Sum all the images in the collection
var sum_image = ee.ImageCollection(mosaic_coll).sum().toInt();


//  Export
Export.image.toDrive({
  image: sum_image, 
  description: 'percountry_LAO_firepolygoncount_2012',
  folder: 'firepolygoncountperpixel',
  fileNamePrefix: 'percountry_LAO_firepolygoncount_2012',
  region: country, //box
  scale: 30.000000000001137, // This is the scale of the original land cover map
  crs: 'EPSG:4326',
  maxPixels: 1e13,
  fileFormat: 'GeoTIFF'
});

However, this task had a runtime of 2 days (and was attempted 2 times, as shown in the task details). I am afraid that this waiting time is too long for me, as I need to repeat this process for the other four countries of mainland SE Asia and also six more times for six more years (currently the script is running only for the year 2012).

My hunch is that the bottleneck is in the Export.toImage task, as it is a client-side function. How do I speed up this process? I already ensured that I used .map() instead of a for-loop. Can anyone suggest modifications to the code? Or should I reach out to Google Earth Engine to request for additional computational power since this may be a heavy geoprocessing task?

Please see below a screenshot of the GEE profiler.

enter image description here

6
  • Do I understand correctly that you basically want to know how many pixels in the fire polygons were affected by fire? So basically, over a period of 1 year, what is the area of those polygons that burned?
    – JonasV
    Nov 15 '21 at 9:20
  • Your issue seems to be that the fire polygons feature collection is way too large. If you have the same data but as a raster, I would use it as a raster. Even just printing the variable polygons is taking ages for me.
    – JonasV
    Nov 15 '21 at 9:58
  • No, that is not what I am trying to achieve. Essentially, I would like a "fire polygon density" map, where each pixel has a value which corresponds to the number of time it intersects with a fire polygon. Nov 16 '21 at 11:10
  • Hmm, does that mean the only alternative for me is to break up my tasks into tiles? Although that would mean that I would have to do a post-hoc merge of all the tiles. I tried setting a shardSize parameter for the export task, but it did not work. Was hoping that GEE would be able to handle such a computation. Nov 16 '21 at 11:11
  • Where does your fire polygon come from? It's pixels converted to polygons right?
    – JonasV
    Nov 16 '21 at 11:16
3

I exported and downloaded as shapefile your polygons vector layer. As it can be observed in following image with QGIS, it has not an id in its attributes table; totaling 127,433 features.

enter image description here

So, I added one and imported to my assets an arbitrary little selection of 27 features, as follows, for exploring where issues could be presents.

enter image description here

I tried out your code with this very low quantity of features and it was impossible to get an useful image for sum_image. So, issues probably are related to firepolygoncountperpixel function.

For this reason, I modified your code and ran it by using the following selection of 27,174 features.

enter image description here

Code looks as follows:

//  Insert MSEA country boundaries
var msea = ee.FeatureCollection('users/jjohanness1992/GADM/mainland_SEA');
var country = msea.filter(ee.Filter.eq('ISO','LAO'));
//print('country metadata',country);

//  Load image
var image = ee.Image('users/jjohanness1992/Fire_and_LCC/clip_lcc_2011_2012_templateLCC').clip(country);
// this is a 30 m raster for the entire mainland SE Asia

//  Create reclassification codes
//  Here I reclassify the land cover map to a raster of 0's
//  This is because I want to preserve this raster grid when performing the polygon count
var oldgroup = ee.List([18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306
]);
var newgroup = ee.List([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
]);

//  Reclassify lcc image to an image of 0 values
var image2 = image.remap(oldgroup,newgroup).select('remapped').rename('b1');

var pol_list = polygons4.toList(polygons4.size());

var coll = ee.ImageCollection(pol_list.map(function (ele){
  
  return image2.add(1)
          .clip(ele);
  
})).sum().toInt();

print(coll);

Map.centerObject(polygons1);
//Map.addLayer(coll);

Export.image.toDrive({
  image: coll, 
  description: 'percountry_LAO_firepolygoncount_2012',
  folder: 'GEE_Folder',
  region: polygons4, //box
  scale: 30, // This is the scale of the original land cover map
  maxPixels: 1e13,
});

After running above code in GEE code editor, exported image (in only 29 minutes with my slow speed Internet service) with Zoom in an area where features represent a very high overlapping, looks as follows. Cursor is pointing out a pixel with a value of 35 in blue area. Red areas represent zero overlapping.

enter image description here

Below image results when it is marked polygons4. It can be corroborated that red areas represent no overlapping.

enter image description here

3
  • Hi @xunilk, your approach and results looks promising. Could you share the polygon assets you used in your script with me? When I open the link to your GEE script, I am unable to run it because the assets are not shared. I am asking this because even after adapting your script to mine, no export task appeared. Thus, I wanted to confirm using your script if an export task can be generated, or if it is something regarding my newly ingested polygon assets (with an integer ID per feature as you suggested). Nov 22 '21 at 15:39
  • 1
    I shared my asset. Try following link: code.earthengine.google.com/8c9ee5f1a5c2beefd6228909c6fd9231
    – xunilk
    Nov 22 '21 at 18:55
  • Thank you @xunilk. I was able to confirm that the link to your new GEE script works. I will try it out on my end again and update this thread accordingly. In addition, could you explain why you converted the ee.FeatureCollection (i.e., polygons4) into a list (i.e., pol_list)? It seemed to speed up the computation, but I found this confusing as I thought an ee.List does not contain coordinate or geometry information like in an ee.FeatureCollection or ee.Feature. Nov 23 '21 at 17:29

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.