6

given is a raster as well as a SpatialPolygonsDataframe. In order to retrieve the highest value of the raster within the area of a polygon, raster::extract can be used. It works fine.

How to get additionally the coordinates of the extracted highest value of the raster within the area of a polygon?

# create raster
r <- raster(ncol=36, nrow=18)
r[] <- runif(ncell(r))
# create SpatialPolygons from GridTopology
grd <- GridTopology(c(-150, -50), c(40, 40), c(8, 3))
Spol <- as(grd, "SpatialPolygons")
# create SpatialPolygonsDataFrame
centroids <- coordinates(Spol)
x <- centroids[,1]
y <- centroids[,2]
SPDF <- SpatialPolygonsDataFrame(Spol, data=data.frame(x=x, y=y, row.names=row.names(Spol)))
# extract max value of raster for each SpatialPolygon
ext <- raster::extract(r, SPDF, fun=max)

*example code is taken from R-documentation

N'ya
  • 347
  • 3
  • 13

2 Answers2

8

You can use the cellnumbers=TRUE argument in extract. With "terra" you can do

library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
v <- vect(system.file("ex/lux.shp", package="terra"))
e <- terra::extract(r, v, xy=TRUE, ID=TRUE)
m <- merge(aggregate(elevation~ID, max, data=e), e)
head(m)
#  ID elevation        x        y
#1  1       547 6.020833 50.17917
#2 10       432 5.987500 49.46250
#3 11       427 6.204167 49.67917
#4 12       413 6.154167 49.77083
#5  2       514 6.087500 49.97917
#6  3       517 5.879167 49.84583

Another way to do the last step

b <- by(e, e$ID, FUN=function(d) d[which.max(d$elevation),])
b <- do.call(rbind, b)
head(b)
#  ID elevation        x        y
#1  1       547 6.020833 50.17917
#2  2       514 6.087500 49.97917
#3  3       517 5.879167 49.84583
#4  4       520 6.112500 49.97917
#5  5       511 5.937500 49.87917
#6  6       403 6.270833 49.83750

And here with your example data and the older packages:

ext <- raster::extract(r, SPDF, cellnumbers=TRUE)
v <- t(sapply(ext, function(i) i[which.max(i[,2]), ] ))

#      cell     value
# [1,]  185 0.9303460
# [2,]  188 0.9821190
# [3,]  154 0.9926290
# [4,]  232 0.8907819
# [5,]  234 0.9998510

To get the coordinates:

xyFromCell(r, v[,1])

#         x   y
# [1,] -135  35
# [2,] -105  35
# [3,]  -85  45
# [4,]  -25  25
# [5,]   -5  25
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
0

A more robust and contemporary solution can be proposed using the terra package and dplyr library. Instead of generating random data, we can use the files from the terra package.

library(terra)   # CRAN v1.7-29
library(dplyr)   # CRAN v1.1.2
library(mapview) # [github::r-spatial/mapview] v2.11.0.9006

# Create a SpatRaster from a file
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

# SpatVector from file
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)

# Extract all values within every polygon
# ID is the polygon unique number
ext <-
  terra::extract(
    r,
    v,
    xy = TRUE,
    ID = TRUE
  )

The resulting dataframe has the polygon ID, coordinates (x, y), and raster value as columns.

>   head(ext)
  ID elevation        x        y
1  1       547 6.020833 50.17917
2  1       485 5.970833 50.17083
3  1       497 5.979167 50.17083
4  1       515 5.987500 50.17083
5  1       515 5.995833 50.17083
6  1       515 6.004167 50.17083

Next, we need to retain only the maximum elevation within each group and transform it into points.

max_p <-
    ext |>
    dplyr::group_by(ID) |>
    dplyr::filter(elevation == max(elevation, na.rm = T)) |>
    terra::vect(geom=c("x", "y"), crs = terra::crs(r))

Finally, we can plot the results interactively with the mapview package.

  mapview::mapview(v, alpha.region = .2) +
    mapview::mapview(max_p, color = "red4", cex = 4, lwd = 2)
atsyplenkov
  • 1,158
  • 13
  • 25