0

I am trying to compute the distance to the nearest waterbodies at a global scale. I am defining waterbodies by thresholding the JRC's Global Surface Water "occurrence" dataset at 80% or greater. This is easy enough to do with fastDistanceTransform(), but the problem is that I need to perform the distance transform at the native scale of the GSW data (30 meters). Ultimately I want to export at a scale of 1 km, but I want the resampling to occur on the distance-transformed image, not the GSW data. This is because coarsening the GSW data can lead to the disappearance of small waterbodies, which greatly alters the final distance transform.

So to be clear, I'd like to:

  1. Threshold GSW at 80% to identify water pixels.
  2. Compute distance transform on thresholded image from (1) at a 30m resolution
  3. Resample (2) at a 1 km resolution for export.

Here's some an example code that fails:

var globe = ee.Geometry.Polygon([-180, 88, 0, 88, 180, 88, 180, -88, 0, -88, -180, -88], null, false)

// Need to perform the distance calculation at native (30m) resolution, then reproject the result. 
var I = JRC.gte(80)
I = I.reproject({crs: 'EPSG:4326',
                 scale: 30})
I = I.fastDistanceTransform({neighborhood:1024}).sqrt().multiply(ee.Image.pixelArea().sqrt())
I = I.reproject({crs: 'EPSG:4326',
                 scale: 1000})

//Map.addLayer(I)
Export.image.toDrive({
  image: I,
  region: globe,
  maxPixels: 1e13,
  folder: 'GEexports',
  description: 'nearestwater_1k',
  scale: 1000
});

You can see that I tried to reproject the image and specify a scale to force the analysis at 30m. I get the error message

Error: Number of pixels requested from Image.reproject exceeds the maximum allowed (2^31).

which is raised by the first reproject (to a scale of 30).

The problem is that if I don't do this reprojection and only specify scale in the Export command, I'm fairly certain GEE will resample the GSW image before thresholding, which is what I'd like to avoid. Nearest distance calculations in this case are highly sensitive to the resolution.

Here's a link to the script.

(I also note that I was able to do the calculation and export at 90 m scale but it produces ~200 files that sum to ~100 GB and I'd like to avoid downloading them, stitching them together into a VRT, then resampling locally if possible.)


I wanted to demonstrate the difference between (1) resampling the GSW image to 1000 meters, thresholding it, and then applying the FDT() and (2) leaving the GSW image in its 30 m (native) resolution, applying the FDT(), and then resampling that output to a 1000 m resolution.

enter image description here

You can see that there is a very large discrepancy between the methods--this is because resampling the GSW dataset to 1000 m basically hides many of the surface water pixels, so the distance transform doesn't "see" many waterbodies of the coarser GSW image. (Top-right resampling was done via gdal on my machine. Bottom-right is GEE output.)

2 Answers 2

3

Exporting at a scale of 1000m means that each output tile has an edge size of 1000m * 1024 pixels = 1024km, but you're using a neighborhood of 1000 pixels in the FDT, so that's a total edge length of 3072km of input to compute each output tile. At the original 30m/pixel, that's an input tile of >10 billion pixels, per output tile.

You probably meant to use a neighborhood of 1024m, not 1024km. Just use a smaller neighborhood.

8
  • I thought exporting at a scale of 1000m meant that each output pixel (not tile?) has a resolution of 1km. If the first reproject works, and therefore forces the distance transfer to operate on a 30m resolution, wouldn't the neighborhood be 1024 pixels x 30 m ~ 30 km? I'm sure I'm just misunderstanding something here...
    – Jon
    Aug 23, 2021 at 14:42
  • The output tile is 1024 pixels, period. And you've specified that the output should be at a resolution of 1000m/pixel, so that tile covers 1024km x 1024km. The neighborhood in the FDT will be measured in output pixels (or more accurately, in the first projection specified after that call, which is the 1km reproject), so that expands the neighborhood out to a total size of 3072 km. To create that output tile, all the required inputs are collected, which, as mentioned, is ~10B pixels. Aug 23, 2021 at 17:18
  • I appreciate the help here, but I'm unable to parse an answer to the question: is there a way to compute the FDT at 30 m scale but export the result at 1000 m scale? I tried reducing the neighborhood, and while I get a successful export, it appears to me that it is computing FDT on the thresholded GSW image at a scale of 1000 m.
    – Jon
    Aug 23, 2021 at 18:14
  • See update for the difference in resampling before/after FDT.
    – Jon
    Aug 23, 2021 at 22:18
  • Lets try this the other way around: how big of a neighborhood, in meters, are you hoping to use for the FDT? Aug 24, 2021 at 19:39
0

Noel provided the real answer in the comments to his answer.

So you just want to move the first reproject to after the FDT: code.earthengine.google.com/225f1ac9bb673c00b57bf677b0c729b8 Remember, reprojections work backwards from where they're called: https://developers.google.com/earth-engine/guides/projections#reprojecting

I can confirm that this export matched my local resampling shown in the upper-right panel of the OP.

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.