CODEWITHSUNDEEP
rrasterpoint
user18018096
user18018096

Reputation:

Finding the peak of a mountain

so I've combined those 2 rasters and made them into one dem raster which contains elevation values:

dem1 = read_stars("srtm_43_06.tif")
dem2 = read_stars("srtm_44_06.tif")
pol = st_read("israel_borders.shp")
dem = st_mosaic(dem1, dem2)
dem = dem[, 5687:6287, 2348:2948]
names(dem) = "elevation"
dem = st_warp(src = dem, crs = 32636, method = "near", cellsize = 90)

Now I need to calculate a point geometry of the peak of the mountain by finding the centroid of the pixel that has the highest elevation in the image, does anyone know what functions I can use?

Upvotes: 1

Views: 238

Answers (2)

Grzegorz Sapijaszko
Grzegorz Sapijaszko

Reputation: 3604

Let's use terra, however similar approach can be applied by raster package as well. For testing purposes we will use raster supplied with terra package

library(terra)
#> terra 1.5.12
f <- system.file("ex/elev.tif", package="terra")
v <- rast(f)
plot(v)

You can check the details of your raster just typing the raster object name and pressing enter, you can check the min and max values with minmax() function form terra:

minmax(v)
#>      elevation
#> [1,]       141
#> [2,]       547

Let's create another raster by copying original one, however checking if the value is the max value of elevation:

w <- v == minmax(v)[2]
plot(w)

Let's create a substitution matrix, and substitute all FALSE with NA and TRUE with 1:

mx <- matrix(c(FALSE, NA, TRUE, 1), ncol = 2, byrow = TRUE)
w <- classify(w, mx)

plot(v)
plot(as.polygons(w), add=TRUE)

Let's find centroids of those polygon(s):

pts <- centroids(as.polygons(w))
plot(pts, add=TRUE)

Let's see our coordinates:

as.data.frame(pts, geom = "WKT")
#>   elevation                   geometry
#> 1         1 POINT (6.020833 50.179167)

Created on 2022-01-29 by the reprex package (v2.0.1)

Upvotes: 1

Robert Hijmans
Robert Hijmans

Reputation: 47156

Building on Grzegorz Sapijaszko's example, here is an alternative path to the top of the mountain.

library(terra)
f <- system.file("ex/elev.tif", package="terra")
x <- rast(f)

If there is a single maximum, you can do

g <- global(x, which.max)
xyFromCell(x, g[,1])
#            x        y
#[1,] 6.020833 50.17917

Now, consider a situation with multiple maxima. I add three more cells with the maximum value.

x[c(1000, 2500, 5000)] <- 547

We can find the four highest peaks with:

g <- global(x, which.max)[[1]]
v <- x[g] |> unlist()
y <- ifel(x == v, v, NA)
p <- as.points(y)
crds(p)
#[1,] 6.020833 50.17917
#[2,] 6.154167 50.10417
#[3,] 5.987500 49.97083
#[4,] 6.237500 49.75417

You should not warp (project with terra) the raster data first because that changes the cell values and potentially the location of the highest peak. You should find the peaks with the original data, but then you can transform the results like this.

pp <- project(p, "EPSG:32636") 
crds(pp)
#            x       y
#[1,] -1411008 5916157
#[2,] -1404896 5904422
#[3,] -1422145 5894509
#[4,] -1413735 5864236

With your files, you could start with something like

ff <- c("srtm_43_06.tif", "srtm_44_06.tif")
v <- vrt(ff)
g <- global(x, which.max)

And then continue as in the examples above.

Upvotes: 2

Related Questions