Inverse distance weighting based on direction

7 views (last 30 days)
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.

Accepted Answer

Wolfgang Schwanghart
Wolfgang Schwanghart on 21 Feb 2022
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)
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 Schwanghart
Wolfgang Schwanghart on 22 Feb 2022
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);
Here's the image output
as well as a detail:
Cheers, Wolfgang

Sign in to comment.

More Answers (0)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by