0

I have a raster (see the plot below) that was created by masking from a bigger raster using a line shapefile. The final aim is to remove spikes and valleys and smooth the raster values along the line. But the problem is I can't figure out how to first sort the values based on the order of occurrence when tracing the path from one end to another. As you can see I can't sort by lat or long as the path doesn't follow a single direction. Obviously, I can't sort by using the cell values (as there are spikes and valleys along the way).

Is there any algorithm/method to number raster cells based on the order of occurrence along a path?

I am mainly using R, so any suggestion using R functions would be great.

Note: I asked a similar question in StackOverflow, but no answers yet. I feel this question is more relevant to this site, so posting a slightly modified version of the question.

enter image description here

Data (in R format)

raster::rasterFromXYZ(a) where a is 

structure(list(x = c(332568.9801, 332568.9801, 332578.9801, 332578.9801, 
332588.9801, 332588.9801, 332588.9801, 332598.9801, 332598.9801, 
332598.9801, 332608.9801, 332608.9801, 332608.9801, 332618.9801, 
332618.9801, 332618.9801, 332628.9801, 332628.9801, 332638.9801, 
332638.9801, 332638.9801, 332648.9801, 332648.9801, 332658.9801, 
332658.9801, 332668.9801, 332668.9801, 332678.9801, 332678.9801, 
332688.9801, 332688.9801, 332688.9801, 332688.9801, 332688.9801, 
332688.9801, 332998.9801, 333008.9801, 333018.9801, 333028.9801, 
333038.9801, 333048.9801, 332688.9801, 332698.9801, 332958.9801, 
332968.9801, 332978.9801, 332988.9801, 332998.9801, 333048.9801, 
333058.9801, 333068.9801, 332698.9801, 332708.9801, 332718.9801, 
332748.9801, 332758.9801, 332768.9801, 332778.9801, 332788.9801, 
332798.9801, 332808.9801, 332818.9801, 332828.9801, 332838.9801, 
332848.9801, 332858.9801, 332868.9801, 332878.9801, 332888.9801, 
332898.9801, 332908.9801, 332918.9801, 332928.9801, 332938.9801, 
332948.9801, 332958.9801, 333068.9801, 332718.9801, 332728.9801, 
332738.9801, 332748.9801, 333068.9801, 333068.9801, 333068.9801, 
333068.9801, 333068.9801, 333068.9801, 333078.9801, 333078.9801, 
333078.9801, 333078.9801, 333078.9801, 333078.9801, 333078.9801, 
333088.9801, 333088.9801, 333088.9801, 333088.9801, 333088.9801, 
333088.9801, 333088.9801, 333078.9801, 333088.9801, 333078.9801, 
333078.9801, 333078.9801, 333078.9801, 333078.9801, 333068.9801, 
333078.9801, 333068.9801, 333068.9801, 333068.9801, 333058.9801, 
333068.9801, 333058.9801, 333058.9801, 333068.9801, 333068.9801, 
333068.9801, 333068.9801, 333078.9801, 333078.9801, 333088.9801, 
333088.9801, 333098.9801, 333108.9801, 333108.9801, 333118.9801, 
333118.9801, 333128.9801, 333138.9801, 333138.9801, 333148.9801, 
333158.9801, 333158.9801, 333168.9801, 333178.9801, 333188.9801, 
333198.9801, 333198.9801, 333208.9801, 333218.9801, 333228.9801, 
333228.9801, 333238.9801, 333248.9801, 333258.9801, 333258.9801, 
333258.9801, 333268.9801, 333268.9801, 333278.9801, 333278.9801, 
333278.9801, 333288.9801, 333288.9801, 333298.9801, 333298.9801, 
333308.9801, 333308.9801, 333318.9801, 333318.9801, 333328.9801, 
333328.9801, 333338.9801, 333338.9801, 333338.9801, 333348.9801, 
333348.9801, 333358.9801), y = c(3061197.4038, 3061187.4038, 
3061187.4038, 3061177.4038, 3061177.4038, 3061167.4038, 3061157.4038, 
3061157.4038, 3061147.4038, 3061137.4038, 3061137.4038, 3061127.4038, 
3061117.4038, 3061117.4038, 3061107.4038, 3061097.4038, 3061097.4038, 
3061087.4038, 3061087.4038, 3061077.4038, 3061067.4038, 3061067.4038, 
3061057.4038, 3061057.4038, 3061047.4038, 3061047.4038, 3061037.4038, 
3061037.4038, 3061027.4038, 3061027.4038, 3061017.4038, 3061007.4038, 
3060997.4038, 3060987.4038, 3060977.4038, 3060977.4038, 3060977.4038, 
3060977.4038, 3060977.4038, 3060977.4038, 3060977.4038, 3060967.4038, 
3060967.4038, 3060967.4038, 3060967.4038, 3060967.4038, 3060967.4038, 
3060967.4038, 3060967.4038, 3060967.4038, 3060967.4038, 3060957.4038, 
3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 
3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 
3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 
3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 
3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 3060957.4038, 
3060947.4038, 3060947.4038, 3060947.4038, 3060947.4038, 3060947.4038, 
3060937.4038, 3060927.4038, 3060917.4038, 3060907.4038, 3060897.4038, 
3060897.4038, 3060887.4038, 3060877.4038, 3060867.4038, 3060857.4038, 
3060847.4038, 3060837.4038, 3060837.4038, 3060827.4038, 3060817.4038, 
3060807.4038, 3060797.4038, 3060787.4038, 3060777.4038, 3060767.4038, 
3060767.4038, 3060757.4038, 3060747.4038, 3060737.4038, 3060727.4038, 
3060717.4038, 3060707.4038, 3060707.4038, 3060697.4038, 3060687.4038, 
3060677.4038, 3060667.4038, 3060667.4038, 3060657.4038, 3060647.4038, 
3060647.4038, 3060637.4038, 3060627.4038, 3060617.4038, 3060617.4038, 
3060607.4038, 3060607.4038, 3060597.4038, 3060597.4038, 3060597.4038, 
3060587.4038, 3060587.4038, 3060577.4038, 3060577.4038, 3060577.4038, 
3060567.4038, 3060567.4038, 3060567.4038, 3060557.4038, 3060557.4038, 
3060557.4038, 3060557.4038, 3060557.4038, 3060547.4038, 3060547.4038, 
3060547.4038, 3060547.4038, 3060537.4038, 3060537.4038, 3060537.4038, 
3060537.4038, 3060527.4038, 3060517.4038, 3060517.4038, 3060507.4038, 
3060507.4038, 3060497.4038, 3060487.4038, 3060487.4038, 3060477.4038, 
3060477.4038, 3060467.4038, 3060467.4038, 3060457.4038, 3060457.4038, 
3060447.4038, 3060447.4038, 3060437.4038, 3060437.4038, 3060427.4038, 
3060417.4038, 3060417.4038, 3060407.4038, 3060407.4038), ele = c(1201.09436035156, 
1201.66271972656, 1201.96130371094, 1201.96130371094, 1201.96130371094, 
1201.96130371094, 1202.08666992188, 1202.18505859375, 1202.23205566406, 
1202.23205566406, 1202.23205566406, 1202.12902832031, 1202.23205566406, 
1202.23205566406, 1202.130859375, 1202.23205566406, 1202.31982421875, 
1202.41235351562, 1202.41235351562, 1202.66650390625, 1203.02331542969, 
1203.02331542969, 1203.02331542969, 1203.02331542969, 1203.02331542969, 
1203.02331542969, 1203.02331542969, 1203.02331542969, 1203.12805175781, 
1203.12805175781, 1203.27319335938, 1203.27319335938, 1203.27319335938, 
1203.357421875, 1203.71594238281, 1204.94604492188, 1204.94604492188, 
1204.94604492188, 1204.94604492188, 1204.94604492188, 1205.32995605469, 
1203.63305664062, 1203.71594238281, 1205.53466796875, 1205.68115234375, 
1204.92346191406, 1204.94604492188, 1204.94604492188, 1205.32995605469, 
1205.35986328125, 1205.35888671875, 1203.73864746094, 1203.63488769531, 
1203.73864746094, 1203.73864746094, 1204.12939453125, 1204.55810546875, 
1204.55810546875, 1205.43127441406, 1204.63513183594, 1203.80871582031, 
1203.80871582031, 1203.80871582031, 1203.80871582031, 1203.80871582031, 
1203.80871582031, 1204.05310058594, 1204.24133300781, 1204.24133300781, 
1204.24133300781, 1204.24133300781, 1204.412109375, 1205.53466796875, 
1205.46630859375, 1205.53466796875, 1205.53466796875, 1205.36047363281, 
1203.73864746094, 1203.73864746094, 1203.73864746094, 1203.73864746094, 
1205.35888671875, 1205.35888671875, 1205.35888671875, 1205.54125976562, 
1206.54333496094, 1206.54333496094, 1206.54333496094, 1206.78332519531, 
1206.56323242188, 1206.41760253906, 1206.41760253906, 1206.39172363281, 
1206.15344238281, 1206.15344238281, 1206.10241699219, 1206.10400390625, 
1206.10241699219, 1206.10241699219, 1206.10241699219, 1206.10241699219, 
1206.43054199219, 1206.43054199219, 1206.87121582031, 1207.04479980469, 
1207.04479980469, 1207.34387207031, 1207.35961914062, 1207.35961914062, 
1207.35961914062, 1207.35961914062, 1207.28234863281, 1207.41857910156, 
1207.58874511719, 1207.58874511719, 1207.58874511719, 1207.58874511719, 
1207.58874511719, 1207.58874511719, 1207.79296875, 1207.8388671875, 
1207.8388671875, 1207.8388671875, 1207.94274902344, 1207.94274902344, 
1207.89782714844, 1207.94274902344, 1207.94274902344, 1208.01794433594, 
1208.43347167969, 1208.53259277344, 1208.56506347656, 1208.56506347656, 
1208.53442382812, 1208.56506347656, 1208.56506347656, 1208.53625488281, 
1208.76940917969, 1208.76940917969, 1208.76940917969, 1208.76940917969, 
1208.76940917969, 1209.5849609375, 1209.5849609375, 1209.5849609375, 
1209.5849609375, 1209.5849609375, 1209.720703125, 1209.52807617188, 
1209.81579589844, 1209.81579589844, 1209.67565917969, 1209.67565917969, 
1209.67565917969, 1209.67565917969, 1209.67565917969, 1209.67565917969, 
1209.67565917969, 1209.67565917969, 1209.67565917969, 1210.79309082031, 
1210.46276855469, 1210.0283203125, 1210.0283203125, 1210.0283203125, 
1210.0283203125, 1210.03039550781, 1210.0283203125, 1210.0283203125, 
1210.0283203125, 1210.0283203125)), class = "data.frame", row.names = c(NA, 
-171L))
3
  • Have you tried rasterToPolygons and then skeletonisation: stackoverflow.com/questions/9595117/…
    – Spacedman
    Commented Jun 14, 2021 at 20:32
  • Thanks for your comment. But I am not sure skeletonisation will solve the issue. I know the coordinates and cell numbers of the raster cells. I just need to find a way to sort the cells based on the order of occurrence along the path (from one end to another).
    – rm167
    Commented Jun 14, 2021 at 22:54
  • It will give you the path, but yes there's other ways based on nearest neighbours. What have you already tried?
    – Spacedman
    Commented Jun 15, 2021 at 10:55

1 Answer 1

0

Here's the nearest-neighbour approach:

First, setup:

library(raster)
source("path.R") # creates "a"
r = rasterFromXYZ(a)

get x,y,V for the non-NA raster points:

pts = rasterToPoints(r)

this many of them

npts = nrow(pts)

get the nearest k=2 neighbours from the X and Y coords:

library(FNN)
nn2 = get.knn(pts[,1:2], k=2)

build a full edge matrix.

nni = nn2$nn.index
edges = rbind(
    cbind(1:npts,nni[,1]),
    cbind(1:npts,nni[,2])
)

Now connect the edges in a graph:

library(igraph)
g = graph_from_edgelist(edges)

Compute the distances of each node from the first node. This assumes the first node is a starting point.

dpath = distances(g, 1)[1,]

Plot the raster, then the points connected in order of the distance from the first node:

plot(r)
lines(pts[order(dpath),])

enter image description here

There's possibly some minor refinement needed at the start and end because those nodes only have one nearest neighbour that you are interested in, so there's maybe a bit more work needed if these points are not easily identifiable (in this example they are the first and last points in the raster so that's easy).

1
  • Thanks for this. I am trying to understand every line. Meanwhile a quick question - If the dpath is the sorted cell numbers in sequential order why does it contain repeated cell numbers (and a 0) at both ends?
    – rm167
    Commented Jun 15, 2021 at 15:33

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.