5

I've done this with PCRaster and ArcGIS in the past but would like to know if there's a way to get it done in Google Earth Engine with the python API?

Can you use other python libraries like pyDEM or pygeoprocessing with google earth engine objects?

0

3 Answers 3

3

EE does not provide a way to compute flow accumulation out of the box, the reason is that a significant amount of iterative operations is required.

You can either use a 90m global version from the above HAND link (reference). But the problem is that dataset is that the flow accumulation is computed only up to sub-catchment limits, then it resets to zero. An alternative is to use a 90m HydroSHEDS FA generated from LDD, see the following post on the development forum and the script with a flow accumulation dataset generated using this method here: https://code.earthengine.google.com/ae9664cd08e21a474f7c93fa6c0d325d.

1

I've experimented a bit with doing this and indeed as Gennaddii says, you'll need a significant amount of iterative operations.

The script below accumulates the values in the image called "data" based on the flow directions given by "Hydrosheds/30DIR". This works, except that it's not really possible to do sufficient iterations to make it useful (in my opinion), but I'll just post it here in case someone else has some ideas on this (https://code.earthengine.google.com/ba976e61fe8cc75062494750d99edb10):

//////////////////////////////////////////
//////////     INPUTS     ////////////////
//////////////////////////////////////////

// Set the number of routing iterations:
var iteration_steps = ee.Number(120);

// Give an image with some data to rout:
var data = ee.Image(1);
// or E.G.: data = ee.Image.pixelArea().divide(1e6);

// Give the bandname from 'data' to rout:
var bandname = ee.String("constant");
// or E.G.: var bandname = ee.String("area");

// Give unit of the data:
var unit = "pixels";
// or E.G.: unit = "km2";

// Give image with flow directions:
var dir = ee.Image("WWF/HydroSHEDS/30DIR").select(["b1"]);

// Set Area-Of-Interest:
var geom = geometry1;

// Set output data type:
var typ = "uint32";

//////////////////////////////////////////
//////////////////////////////////////////
//////////////////////////////////////////

// Rename the band.
data = data.select([bandname],["rout"]);

// Add a band to store the summation of the routs.
data = data.addBands(data.select(["rout"],["summed"]));

// Set property indicating the current iteration step.
data = data.set("iter_idx", 0);

// Create a feature-collection that describes to what directions the numbers in 'dir' refer to (e.g. 128 is NE).
var col = ee.FeatureCollection(ee.List([
          ee.Feature(null, {"weight": ee.List([[0,0,0],[1,0,0],[0,0,0]]),
                            "direction": "E",
                            "dir_id": 1,}),
          ee.Feature(null, {"weight": ee.List([[1,0,0],[0,0,0],[0,0,0]]),
                            "direction": "SE",
                            "dir_id": 2,}),
          ee.Feature(null, {"weight": ee.List([[0,1,0],[0,0,0],[0,0,0]]),
                            "direction": "S",
                            "dir_id": 4,}),
          ee.Feature(null, {"weight": ee.List([[0,0,1],[0,0,0],[0,0,0]]),
                            "direction": "SW",
                            "dir_id": 8,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,1],[0,0,0]]),
                            "direction": "W",
                            "dir_id": 16,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,0],[0,0,1]]),
                            "direction": "NW",
                            "dir_id": 32,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,0],[0,1,0]]),
                            "direction": "N",
                            "dir_id": 64,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,0],[1,0,0]]),
                            "direction": "NE",
                            "dir_id": 128,})]));

// The iteration algorithm.
var iterate_route = function(iter_step, data){

    // The function to do one routing iteration for one direction.
    var route = function(ft){

        // Create the kernel for the current flow direction.
        var kernel = ee.Kernel.fixed({width: 3, 
                                      height: 3, 
                                      weights: ft.get("weight")
                                     });

        // Get the number corresponding to the flow direction map.
        var dir_id = ee.Number(ft.get("dir_id"));

        // Mask irrelevent pixels.
        var routed = ee.Image(data).select("rout").updateMask(dir.eq(dir_id));

        // Move all the pixels one step in the direction currently under consideration.
        var routed_dir = routed.reduceNeighborhood({reducer: ee.Reducer.sum(), 
                                                    kernel: kernel,
                                                    skipMasked: false
                                                   });

        return routed_dir;
        };

    // Loop over the eight flow directions and sum all the routed pixels together.
    var step = ee.ImageCollection(col.map(route)).sum().rename("rout")
                                                       .reproject(dir.projection())
                                                       .set("iter_idx", iter_step)
                                                       .clip(geom);

    // Sum the newest routed pixels with previous ones.
    var summed = step.select("rout").add(ee.Image(data).select("summed")).rename("summed");

    // Add the 'rout' and 'summed' bands together in one image.
    var data_next_step = step.addBands(summed);

    return data_next_step;
    };

// Create dictionary to cast output to the chosen datatype.
var cast = ee.Dictionary({"rout" : typ,
                          "summed": typ});

// Create the list with iteration step numbers.
var steps = ee.List.sequence(1, iteration_steps);

// Do the actual iterations.
var full = ee.Image(steps.iterate(iterate_route, data)).cast(cast);

// Add Layer to Map
Map.addLayer(full.select("summed"), {min: 0, max: 1000}, "summed");
Map.centerObject(full.select("summed"))
0
0

In case anyone stumbles across this in the future. MERIT Hydro dataset is available on Earth Engine now.

It includes the "upa" band which holds the upstream drainage area of a pixel and has 3 arc seconds resolution.

It is newly added in October 2020. https://medium.com/google-earth/meet-the-newest-additions-to-the-earth-engine-data-catalog-2380a22d63a9

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