Cannot get xy coordinates of a raster object imported with terra #29

Robinlovelace · 2020-04-25T17:41:46Z

Context: I have discovered that raster's extract() function can exactly reproduce the elevation profiles generated by proprietary software in functionality to calculate the gradient, slope, of linestrings. I plan to implement this, as per ropensci/stplanr#392

The problem with the current implementation: it takes an hour to run single threaded on a regional example and could take days on national road datasets. I'm hoping terra could be the tool for the job, reasonable hope?

Reproducible example that I hope supports development efforts (I know this is work in progress):

# Aim: get gradients at xy locations with bilinear interpolation
remotes::install_github("rspatial/terra")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'terra' from a github remote, the SHA1 (a36cf92e) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(terra)
#> This is version 0.5.12 of the "terra" package, for evaluation only
download.file("https://github.com/geocompr/d/releases/download/1/r1.zip", "r1.zip")
unzip("r1.zip")
list.files()
#> [1] "r1"                     "r1.zip"                 "reprex_reprex.R"       
#> [4] "reprex_reprex.spin.R"   "reprex_reprex.spin.Rmd"
list.files("r1")
#> [1] "dblbnd.adf"   "hdr.adf"      "sta.adf"      "w001001.adf"  "w001001x.adf"
dem = terra::rast("r1/hdr.adf")
dem_matrix = terra::xyFromCell(dem, row = 1:terra::nrow(dem), col = 1:terra::ncol(dem))
#> Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'xyFromCell' for signature '"SpatRaster", "missing"'

^{Created on 2020-04-25 by the reprex package (v0.3.0)}

The text was updated successfully, but these errors were encountered:

Robinlovelace · 2020-04-25T17:53:12Z

Another reproducible example that reads-in the vector dataset (terra's read functionality seems very efficient btw). Note: I need estimated elevations at each vertex of the linestring, but was also having issues with the geom() function.

# Aim: get gradients at xy locations with bilinear interpolation
remotes::install_github("rspatial/terra")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'terra' from a github remote, the SHA1 (a36cf92e) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(terra)
#> This is version 0.5.12 of the "terra" package, for evaluation only
download.file("https://github.com/geocompr/d/releases/download/1/r1.zip", "r1.zip")
unzip("r1.zip")
dem = terra::rast("r1/hdr.adf")
download.file("https://github.com/geocompr/d/releases/download/1/example_route_2d.gpkg", "example_route_2d.gpkg", mode="wb")
v = vect("example_route_2d.gpkg")
plot(v)

# terra::extract(dem, v, method = "bilinear") # fails

^{Created on 2020-04-25 by the reprex package (v0.3.0)}

rhijmans · 2020-04-25T21:41:43Z

In the first chunk you probably want dem_matrix = terra::xyFromCell(dem, 1:ncell(dem))
geom worked fine for me. I added an "as.points" method for SpatVector

#dir.create("test", FALSE)
setwd("test")
#download.file("https://github.com/geocompr/d/releases/download/1/r1.zip", "r1.zip")
#download.file("https://github.com/geocompr/d/releases/download/1/example_route_2d.gpkg", "example_route_2d.gpkg", mode="wb")
#unzip("r1.zip")

library(terra)
dem <- terra::rast("r1/hdr.adf")

dem_matrix = terra::xyFromCell(dem, 1:ncell(dem))

v <- vect("example_route_2d.gpkg")

## You may not see the messages below. They can be ignored, but are very annoying
## They are related to garbage collection, but I do not know how to get rid of them
#Error in x$.self$finalize() : attempt to apply non-function
#Error in (function (x)  : attempt to apply non-function
#Error in (function (x)  : attempt to apply non-function

g <- geom(v)

head(g)
#  id part        x         y hole
#1  1    1 -92539.1 -106077.5    0
#2  1    1 -92529.4 -106078.7    0
#3  1    1 -92519.7 -106079.9    0
#4  1    1 -92510.0 -106081.1    0
#5  1    1 -92500.3 -106082.4    0
#6  1    1 -92490.6 -106083.6    0

xy <- g[, c("x", "y")]
head(extract(dem, xy))
#        hdr
#[1,] 90.943
#[2,] 90.462
#[3,] 90.283
#[4,] 90.246
#[5,] 90.200
#[6,] 90.144

# new "as.points" for vector 
p <- as.points(v)

"p" is a multi-point object
extract(dem, p)
#[[1]]
#[[1]][[1]]
# [1] 90.943 90.462 90.283 90.246 90.200 90.144 90.087 90.031 90.000 90.000
#[11] 88.687 87.367

extract(dem, p, method="bilinear")
#[[1]]
#[[1]][[1]]
 #[1] 90.92585 90.46447 90.28246 90.23851 90.18833 90.13170 90.07147 90.02918
 #[9] 89.92668 89.30816 88.25591 87.17434

Robinlovelace · 2020-04-26T08:12:47Z

Thanks Robert. Another issue is that I am seeing this message frequently when using (and at the moment when installing) terra:

But that is probably a separate issue and by the sounds of it you have solved my specific issue. Testing it now...

Robinlovelace · 2020-04-26T08:19:23Z

It's working 🎉

# Aim: get gradients at xy locations with bilinear interpolation
# See https://github.com/rspatial/terra/issues/29
remotes::install_github("rspatial/terra")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'terra' from a github remote, the SHA1 (35f7d051) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(terra)
#> This is version 0.6.1 of the "terra" package, for evaluation only
download.file("https://github.com/geocompr/d/releases/download/1/r1.zip", "r1.zip")
unzip("r1.zip")
dem = terra::rast("r1/hdr.adf")
download.file("https://github.com/geocompr/d/releases/download/1/example_route_2d.gpkg", "example_route_2d.gpkg")
v = vect("example_route_2d.gpkg")
plot(v)

g = geom(v)
head(g)
#>   id part        x         y hole
#> 1  1    1 -92539.1 -106077.5    0
#> 2  1    1 -92529.4 -106078.7    0
#> 3  1    1 -92519.7 -106079.9    0
#> 4  1    1 -92510.0 -106081.1    0
#> 5  1    1 -92500.3 -106082.4    0
#> 6  1    1 -92490.6 -106083.6    0
z = terra::extract(dem, g[, c("x", "y")], method = "bilinear") # fails
z
#>          hdr
#>  [1,] 90.943
#>  [2,] 90.462
#>  [3,] 90.283
#>  [4,] 90.246
#>  [5,] 90.200
#>  [6,] 90.144
#>  [7,] 90.087
#>  [8,] 90.031
#>  [9,] 90.000
#> [10,] 90.000
#> [11,] 88.687
#> [12,] 87.367

dem_matrix = xyFromCell(dem, 1:ncell(dem))
head(dem_matrix)
#>            x      y
#> [1,] -113450 -79570
#> [2,] -113440 -79570
#> [3,] -113430 -79570
#> [4,] -113420 -79570
#> [5,] -113410 -79570
#> [6,] -113400 -79570

^{Created on 2020-04-26 by the reprex package (v0.3.0)}

Nowosad · 2020-04-26T08:23:53Z

@Robinlovelace idea - create some benchmarks

Robinlovelace · 2020-04-26T08:30:09Z

Equivalent code in raster with timings

# Aim: get gradients at xy locations with bilinear interpolation
library(raster)
#> Loading required package: sp
download.file("https://github.com/geocompr/d/releases/download/1/r1.zip", "r1.zip")
unzip("r1.zip")
tictoc::tic()
dem = raster("r1/hdr.adf")
tictoc::toc()
#> 0.072 sec elapsed
download.file("https://github.com/geocompr/d/releases/download/1/example_route_2d.gpkg", "example_route_2d.gpkg")
tictoc::tic()
v = rgdal::readOGR("example_route_2d.gpkg")
#> Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
#> dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in CRS
#> definition: +proj=tmerc +lat_0=39.6682583333333 +lon_0=-8.13310833333333 +k=1
#> +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> OGR data source with driver: GPKG 
#> Source: "/tmp/RtmpUnnj6J/reprex58e81cf15c2c/example_route_2d.gpkg", layer: "example_route_2d"
#> with 1 features
#> It has 0 fields
plot(v)

g = geom(v)
head(g)
#>      object part cump        x         y
#> [1,]      1    1    1 -92539.1 -106077.5
#> [2,]      1    1    1 -92529.4 -106078.7
#> [3,]      1    1    1 -92519.7 -106079.9
#> [4,]      1    1    1 -92510.0 -106081.1
#> [5,]      1    1    1 -92500.3 -106082.4
#> [6,]      1    1    1 -92490.6 -106083.6
tictoc::toc()
#> 0.192 sec elapsed
tictoc::tic()
z = extract(dem, g[, c("x", "y")], method = "bilinear") # fails
tictoc::toc()
#> 0.1 sec elapsed
tictoc::tic()
dem_matrix = xyFromCell(dem, 1:ncell(dem))
tictoc::toc()
#> 0.258 sec elapsed
head(dem_matrix)
#>            x      y
#> [1,] -113450 -79570
#> [2,] -113440 -79570
#> [3,] -113430 -79570
#> [4,] -113420 -79570
#> [5,] -113410 -79570
#> [6,] -113400 -79570

^{Created on 2020-04-26 by the reprex package (v0.3.0)}

Robinlovelace · 2020-04-26T08:31:18Z

Profile output for raster:

Robinlovelace · 2020-04-26T08:34:47Z

Result: terra seems faster on my set-up.

# Aim: get gradients at xy locations with bilinear interpolation
    # See https://github.com/rspatial/terra/issues/29
    remotes::install_github("rspatial/terra")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'terra' from a github remote, the SHA1 (35f7d051) has not changed since last install.
#>   Use `force = TRUE` to force installation
    library(terra)
#> This is version 0.6.1 of the "terra" package, for evaluation only
    download.file("https://github.com/geocompr/d/releases/download/1/r1.zip", "r1.zip")
    unzip("r1.zip")
    tictoc::tic()
    dem = terra::rast("r1/hdr.adf")
    tictoc::toc()
#> 0.003 sec elapsed
    download.file("https://github.com/geocompr/d/releases/download/1/example_route_2d.gpkg", "example_route_2d.gpkg")
    tictoc::tic()
    v = vect("example_route_2d.gpkg")
    plot(v)

    g = geom(v)
    head(g)
#>   id part        x         y hole
#> 1  1    1 -92539.1 -106077.5    0
#> 2  1    1 -92529.4 -106078.7    0
#> 3  1    1 -92519.7 -106079.9    0
#> 4  1    1 -92510.0 -106081.1    0
#> 5  1    1 -92500.3 -106082.4    0
#> 6  1    1 -92490.6 -106083.6    0
    tictoc::toc()
#> 0.062 sec elapsed
    tictoc::tic()
    z = terra::extract(dem, g[, c("x", "y")], method = "bilinear") # fails
    tictoc::toc()
#> 0.003 sec elapsed
    tictoc::tic()
    dem_matrix = xyFromCell(dem, 1:ncell(dem))
    tictoc::toc()
#> 0.99 sec elapsed
    head(dem_matrix)
#>            x      y
#> [1,] -113450 -79570
#> [2,] -113440 -79570
#> [3,] -113430 -79570
#> [4,] -113420 -79570
#> [5,] -113410 -79570
#> [6,] -113400 -79570

^{Created on 2020-04-26 by the reprex package (v0.3.0)}

Robinlovelace · 2020-04-26T08:36:06Z

Heads-up @rhijmans xyFromCell() seems slow here FYI:

    tictoc::tic()
    dem_matrix = xyFromCell(dem, 1:ncell(dem))
    tictoc::toc()
#> 0.99 sec elapsed

This was referenced Apr 25, 2020

extent inconsistency #31

Closed

"Error in (function (x) : attempt to apply non-function" #30

Closed

Robinlovelace closed this as completed Apr 26, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Cannot get xy coordinates of a raster object imported with terra #29

Cannot get xy coordinates of a raster object imported with terra #29

Robinlovelace commented Apr 25, 2020

Robinlovelace commented Apr 25, 2020 •

edited by rhijmans

rhijmans commented Apr 25, 2020 •

edited

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Nowosad commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Cannot get xy coordinates of a raster object imported with terra #29

Cannot get xy coordinates of a raster object imported with terra #29

Comments

Robinlovelace commented Apr 25, 2020

Robinlovelace commented Apr 25, 2020 • edited by rhijmans

rhijmans commented Apr 25, 2020 • edited

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Nowosad commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 26, 2020

Robinlovelace commented Apr 25, 2020 •

edited by rhijmans

rhijmans commented Apr 25, 2020 •

edited