2

I am working in GEE and I am trying to extract areas of clusters for some specific coordinates.

I have a clustered map of the area of all the clusters (using object-based clustering) and a list of coordinates. I am trying to extract areas of those clusters only whose coordinates are available. The coordinates can lie anywhere in the cluster but that cluster's whole area should be extracted.

I have tried this approach but it isn't giving me the area of the whole cluster. What I am getting are the values like: 11.8, 8.7, etc. But the actual values are in 100's.

area_map = ee.Image.pixelArea().divide(1e6).addBands(clusters).reduceConnectedComponents(ee.Reducer.sum(), "clusters", 256) #Per-cluster area calculation
coordinates_features = ee.FeatureCollection(coordinates.tolist())
area_list = area_map.reduceRegions(collection= coordinates_features, reducer= ee.Reducer.first().setOutputs(['area']), scale=30)
print('area', area_list)

Can someone please guide, how do we achieve this?

2
  • Please, edit your question and put your own code. It is not clear which is the relation between the accepted answer with your question.
    – xunilk
    May 30, 2023 at 21:35
  • @xunilk I have added the code chunk now.
    – th145
    May 31, 2023 at 5:08

1 Answer 1

2
+50

I'm not sure I understand what you want to do completely. It would have been helpful if you provided a complete, runnable script.

Here's my interpretation of what you're asking, correct me if I got this wrong. You have a cluster image, with a categorical image band, and a set of points. You want to find the area of the cluster each of these points fall on.

reduceConnectedComponents() have a max size, and when you make that value too large, you'll run out of memory. So you cannot find the complete area of a cluster with this approach.

Here's an idea:

  • Find the clusters you're interested in by sampling the cluster image at your points
  • Mask out the part of the cluster image with other cluster ids
  • reduceRegions() on the pixel area and the masked cluster image, use a sum reducer grouped by the cluster id.

In principle, I suppose you could skip the first two steps too, you'd get the area for every cluster id. But if you have a lot of them, maybe it's better to skip them.

Here it is in JavaScript - it should be trivial to port to Python:

var clusterIds = clusters
  .sampleRegions({collection: points, scale: 10})
  .aggregate_array('cluster')
  .distinct()
var mask = ee.Image(ee.Array(clusterIds)).eq(clusters)
  .arrayReduce(ee.Reducer.max(), [0])
  .arrayFlatten([['cluster']])
var maskedClusters = clusters.updateMask(mask)
var areas = ee.List(
  ee.Image.pixelArea()
    .addBands(maskedClusters)
    .reduceRegion({
      reducer: ee.Reducer.sum().group(1), 
      geometry: region, 
      scale: 10,
      maxPixels: 1e13
    })
    .get('groups')
)   

https://code.earthengine.google.com/48cd99979238a57222b63def0351d9b2

7
  • Thank you for your answer. Sorry for the confusion, actually I am getting areas of the clusters (in km^2) but they are incorrect. I am using object-based clustering. The areas that I am getting are like: 11.8, 8.7, etc. But the actual areas of the clusters are in 100's. I ran your code but your code isn't giving me the areas (in numbers) at all.
    – th145
    Jun 2, 2023 at 14:02
  • I've updated my answer. I had made a mistake with the area calculation. I cannot help you more than this unless you providing more details.. Jun 5, 2023 at 8:53
  • Thank you for you help Daniel. I am using this repo by Noel Gorelick: code.earthengine.google.com/?accept_repo=users/gorelick/EE102 Could you please explain what I am doing wrong here?
    – th145
    Jun 5, 2023 at 11:46
  • the changes aren't getting saved in this repo. I added your code chunk to get the area. But I am getting this error: List (Error) Array: No numbers in 'values', must provide a type. mask: Layer error: Array: No numbers in 'values', must provide a type. maskedClusters: Layer error: Array: No numbers in 'values', must provide a type. Here are the coordinates: [-121.7142,38.9714], [-121.8845,38.9735], [-121.8832,38.9297] and [-121.7211,38.9308]
    – th145
    Jun 5, 2023 at 11:52
  • sorry for bothering you. Could you please tell, how can we visualize the "areas" it on map?
    – th145
    Jun 7, 2023 at 9:08

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.