Inverse distance weighting based on direction

This problem is a bit specific, so I'll try to be as clear as I can. I have 1554 points along the coast of Greenland where material is transported into the fjords and sea surrounding the coast. I have a grid, spanning the fjords and sea, consisting of 264192 points, where I would like to know the amount of material transported to each of these. I have made an inverse distance weighting model, where I have made a 1554x264192 matrix containing the euclidean distance from each outlet to each point, and then calculated the weight based on the inverse distance 1/d. This gives me a flat model, only taking the distance from each outlet into consideration. However, a realistic model would only take the outlets coming from upstream into consideration.

I know the bathymetry of the fjords and sea, and can calculate the gradient for each point. The outlets that are further out the sea (with a bigger bathymetry), should not be taken into consideration at any points with a lower bathymetry (further inland). I was thinking that I could implement an if loop, where every outlet at a point downstream should not be taken into consideration, however I'm finding it hard to implement

This is very vague, so I don't expect much help, but I appreciate any. I use Wolfgang Schwangharts topotoolbox, but the dataset is so large that I cannot attach it, unfortunately.

An image of the flat model. The black dots are the outputs from where material is transported, and every point receives material based on the inverse distance.

Thanks in advance.

Answer 1

[Wolfgang Schwanghart]

Hi Christian,

interesting question. First of all, if you follow the IDW-approach, you might save a lot of computing time using the function knnsearch (or rangesearch). This would enable restricting for each bathymetric point the number of outlets.

However, for sediment transported processes, IDW is likely not a good model. Rather, there is most likely a preferred transport direction, and that direction might follow the downward bathymetric gradient (as long as other process such as wave action, shore-parallel currents etc. are not too important).

To model this kind of sediment transport, you might again use TopoToolbox. This would require you to create a GRIDobj (e.g. BATH) from the bathymetric points (possibly by merging it also with shoreline pixels from the topographic DEM). As they are already gridded, this won't be a big issue (otherwise, check the function interp2GRIDobj). Once you have the grid, you can calculate a FLOWobj. I suggest to use multiple flow directions (FLOWobj(fillsinks(BATH),'multi')) which would account for the diffusiveness of fine sediment along their transport path. Then use the locations of outlets (with their coordinates x and y) and determine their linear indices in the BATH grid (IX = coord2ind(BATH,x,y)). Make sure that all BATH.Z(IX) values are non-nan. Otherwise, you will need to move them to the nearest valid pixel in BATH.

Now you can route sediment downslope entering the sea at the outlets. For this, create a new grid W (e.g. with W = GRIDobj(BATH)) and set values in W.Z(IX) = 1; You now have a grid with zeros and ones at outlets that you can use in a weighted flow accumulation approach.

A = flowacc(FD,W)
imagesc(A)

The resulting image shows modelled sediment transport paths given that the amount of sediment entering the sea at each outlet point is equal. Of course, you can also modify the approach to give heigher weights to these outlets with larger contributing area and likely higher sediment input.

Hope this helps. Best regards,

Wolfgang

[Christian Mathiesen]

Hi Wolfgang

Thank you for your response, I appreciate it.

I should've mentioned that the DEM also contains the topographic values of inland Greenland, and not only the bathymetry, which is how I was able to locate the outlets using STREAMobj and klargestconncomps. However, when I calculate A = flowacc(FD), it appears to go inland as well, as seen on the figure below. Is there a way to circumvent this, as I am only interested in the sediment transport in the sea?

DEM = GRIDobj('DEM.tif'); 
FD  = FLOWobj(fillsinks(DEM)); % I am using single flow until I can use a machine with more RAM
IX = coord2ind(DEM,x_outlet,y_outlet);

W = GRIDobj('DEM.tif');
W.Z(IX) = 1
A = flowacc(FD)
figure
imagesc(A)
hold on
plot(x_outlet,y_outlet,'*k','MarkerSize',5)
hold off
colorbar

Best,

Christian

[Wolfgang Schwanghart]

Hi Christian,

there is a small issue in your code creating the GRIDobj W. And yes, FLOWobj with MFD unfortunately requires quite some RAM.

DEM = GRIDobj('DEM.tif'); 
FD  = FLOWobj(fillsinks(DEM)); % I am using single flow until I can use a machine with more RAM
IX = coord2ind(DEM,x_outlet,y_outlet);
% The following line creates a GRIDobj of zeros.
W = GRIDobj(DEM);
W.Z(IX) = 1;
A = flowacc(FD);
figure
imagesc(A)
hold on
plot(x_outlet,y_outlet,'*k','MarkerSize',5)
hold off
colorbar

Like this, no water should be seen within the terrestrial part of the DEM.

Cheers, Wolfgang

[Christian Mathiesen]

Hi Wolfgang,

That indeed solved it! Now I just have to distribute proper weights to specific outlets with higher/lower sediment input. Thanks a lot for your help.

Best,

Christian.

[Wolfgang Schwanghart]

In terrestrial landscapes, you will often find that Q (discharge) and potentially Qs (sediment discharge) are proportional to sqrt(Upstream Area). I am not sure whether there is such a relation for glaciated landscapes, but it might be a good value to start with.

Cheers, Wolfgang

[Christian Mathiesen]

Just a quick follow-up question.

Will a FLOWobj with MFD distribute flow to the entire fjord, bounded inside by the black lines (drawing the coast)? I expect that there would be distribution of sediment to the whole fjord, and not only along the flowline below. If not, is there another function that can draw additional flow lines?

Regards,

Christian.

[Wolfgang Schwanghart]

That depends on the Fjord's bathymetry. If it's steep and V-shaped (which it is likely not), then MFD will result in very concentrated flow similar to the single (D8) flow. But if you want to have it very dispersive (and more to look like an ice flow model), then you'll need another approach which doesn't (heavily) rely on the bathymetry.

Cheers, Wolfgang

[Wolfgang Schwanghart]

The following code shows how to calculate sediment discharge with directions based on a distance transform.

DEM = readexample('taiwan');
DEM = inpaintnans(DEM);
DEM = pad(DEM,200,nan(1,class(DEM.Z)));

% On land
FD  = FLOWobj(DEM);
S   = STREAMobj(FD,'minarea',1000);
A   = flowacc(FD);
% outlets
IXO = streampoi(S,'out','ix');

% The sea
D   = distance(DEM,erode(~isnan(DEM),ones(3)));
D   = clip(D,D~=0);
D   = 1./D;

FDS  = FLOWobj(fillsinks(D),'multi');
W   = GRIDobj(D);
W.Z(IXO) = A.Z(IXO);
AS  = flowacc(FDS,W);

imageschs(DEM,log(AS+A),'colormap','flowcolor')

Here's the image output

as well as a detail:

Cheers, Wolfgang