search

ropensci / slopes

Package to calculate slopes of roads, rivers and trajectories
https://docs.ropensci.org/slopes/
GNU General Public License v3.0
70 stars 6 forks source link

slope_raster() : Error in weighted.mean.default(abs(g1), d, na.rm = TRUE) : 'x' and 'w' must have the same length #12

Closed temospena closed 3 years ago

temospena commented 3 years ago

Hi,

I'm having a hard time to figure out what am I doing wrong. I have a road network and a DEM. They have the same crs (they overlap), but when using the function slope_raster it retrieves this error.

I'm using a large dataset, available here And a DEM from NASA, with 27m res, downloaded with QGIS' SRTM-downloader plugin, and clipped.

#### Error

#packages
library(sf)
library(raster)
library(slopes)

#import datasets
network = st_read("https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true")
#> Reading layer `RedeViariaLisboa_dadosabertos' from data source `https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true' using driver `GPKG'
#> Simple feature collection with 22098 features and 23 fields
#> geometry type:  MULTILINESTRING
#> dimension:      XYZ
#> bbox:           xmin: -9.229783 ymin: 38.69186 xmax: -9.090931 ymax: 38.79573
#> z_range:        zmin: 0 zmax: 213.8635
#> geographic CRS: WGS 84

network = st_zm(network, drop = T) #make sure it has no Z values stored

dem = raster("https://github.com/U-Shift/Declives-RedeViaria/blob/main/raster/LisboaNASA_clip.tif?raw=true")
dem
#> class      : RasterLayer 
#> dimensions : 430, 729, 313470  (nrow, ncol, ncell)
#> resolution : 0.0002777778, 0.0002777778  (x, y)
#> extent     : -9.27875, -9.07625, 38.68903, 38.80847  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#> source     : https://github.com/U-Shift/Declives-RedeViaria/blob/main/raster/LisboaNASA_clip.tif?raw=true 
#> names      : LisboaNASA_clip.tif.raw.true

#do they overlap?
raster::plot(dem)
plot(sf::st_geometry(network), add = TRUE)


#why this error?
network$slopeNASA = slope_raster(network, e = dem)
#> [1] FALSE
#> Warning in max(d, na.rm = TRUE): no non-missing arguments to max; returning -Inf
#> Warning in max(d, na.rm = TRUE): no non-missing arguments to max; returning -Inf
#> Error in weighted.mean.default(abs(g1), d, na.rm = TRUE): 'x' and 'w' must have the same length

Created on 2020-10-27 by the reprex package (v0.3.0)

Session info ``` r devtools::session_info() #> - Session info --------------------------------------------------------------- #> setting value #> version R version 3.6.2 (2019-12-12) #> os Windows 7 x64 SP 1 #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate Portuguese_Portugal.1252 #> ctype Portuguese_Portugal.1252 #> tz Europe/London #> date 2020-10-27 #> #> - Packages ------------------------------------------------------------------- #> package * version date lib source #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.1) #> backports 1.1.7 2020-05-13 [1] CRAN (R 3.6.3) #> callr 3.4.3 2020-03-28 [1] CRAN (R 3.6.3) #> class 7.3-15 2019-01-01 [1] CRAN (R 3.6.2) #> classInt 0.4-3 2020-04-07 [1] CRAN (R 3.6.3) #> cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.3) #> codetools 0.2-16 2018-12-24 [1] CRAN (R 3.6.2) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.1) #> curl 4.3 2019-12-02 [1] CRAN (R 3.6.2) #> DBI 1.1.0 2019-12-15 [1] CRAN (R 3.6.2) #> desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.1) #> devtools 2.3.0 2020-04-10 [1] CRAN (R 3.6.3) #> digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.3) #> dplyr 1.0.0 2020-05-29 [1] CRAN (R 3.6.3) #> e1071 1.7-3 2019-11-26 [1] CRAN (R 3.6.2) #> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 3.6.3) #> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.1) #> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.2) #> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.1) #> generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.1) #> geodist 0.0.6 2020-10-15 [1] CRAN (R 3.6.3) #> glue 1.4.2 2020-08-27 [1] CRAN (R 3.6.3) #> highr 0.8 2019-03-20 [1] CRAN (R 3.6.1) #> htmltools 0.5.0 2020-06-16 [1] CRAN (R 3.6.3) #> httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.1) #> KernSmooth 2.23-16 2019-10-15 [1] CRAN (R 3.6.2) #> knitr 1.29 2020-06-23 [1] CRAN (R 3.6.3) #> lattice 0.20-38 2018-11-04 [1] CRAN (R 3.6.2) #> lifecycle 0.2.0 2020-03-06 [1] CRAN (R 3.6.3) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.1) #> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.2) #> mime 0.9 2020-02-04 [1] CRAN (R 3.6.2) #> pbapply 1.4-2 2019-08-31 [1] CRAN (R 3.6.1) #> pillar 1.4.4 2020-05-05 [1] CRAN (R 3.6.3) #> pkgbuild 1.0.8 2020-05-07 [1] CRAN (R 3.6.3) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.2) #> pkgload 1.1.0 2020-05-29 [1] CRAN (R 3.6.3) #> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.2) #> processx 3.4.2 2020-02-09 [1] CRAN (R 3.6.3) #> ps 1.3.3 2020-05-08 [1] CRAN (R 3.6.3) #> purrr 0.3.4 2020-04-17 [1] CRAN (R 3.6.3) #> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.2) #> raster * 3.1-5 2020-04-19 [1] CRAN (R 3.6.3) #> Rcpp 1.0.5 2020-07-06 [1] CRAN (R 3.6.3) #> remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.3) #> rgdal 1.4-8 2019-11-27 [1] CRAN (R 3.6.2) #> rlang 0.4.7 2020-07-09 [1] CRAN (R 3.6.3) #> rmarkdown 2.3.8 2020-09-19 [1] Github (rstudio/rmarkdown@fbd0b1f) #> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.1) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.2) #> sf * 0.9-6 2020-09-13 [1] CRAN (R 3.6.3) #> slopes * 0.0.0.9000 2020-05-10 [1] Github (itsleeds/slopes@cdff333) #> sp * 1.4-1 2020-02-28 [1] CRAN (R 3.6.3) #> stringi 1.4.6 2020-02-17 [1] CRAN (R 3.6.2) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.1) #> testthat 2.3.2 2020-03-02 [1] CRAN (R 3.6.3) #> tibble 3.0.1 2020-04-20 [1] CRAN (R 3.6.3) #> tidyselect 1.1.0 2020-05-11 [1] CRAN (R 3.6.3) #> units 0.6-6 2020-03-16 [1] CRAN (R 3.6.3) #> usethis 1.6.1 2020-04-29 [1] CRAN (R 3.6.3) #> vctrs 0.3.1 2020-06-05 [1] CRAN (R 3.6.3) #> withr 2.2.0 2020-04-20 [1] CRAN (R 3.6.3) #> xfun 0.17 2020-09-09 [1] CRAN (R 3.6.3) #> xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.1) #> yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.3) #> #> [1] C:/Program Files/R/R-3.6.2/library ```
Robinlovelace commented 3 years ago

Many thanks for the report @temospena, taking a look now...

Robinlovelace commented 3 years ago

Update: I can reproduce the issue. Looking into it, as documented here: https://github.com/ITSLeeds/slopes/blob/master/vignettes/issues.Rmd

Robinlovelace commented 3 years ago

Think I found the cause of the issue: LINESTRING vs MULTILINESTRING geometries. Converting to the former with st_cast() fixed it. I plan to add a check at the beginning of function that will tell people to convert if it's the wrong geometry type.

Tests:

network_subset = network[1:9, ]
network$slopesNASA = slope_raster(network_subset, e = dem) # also fails with the same error
dem_subset = crop(dem, network)
network_subset$slopesNASA = slope_raster(network_subset, e = dem_subset) # also fails with the same error
unique(st_geometry_type(network_subset))
unique(st_geometry_type(lisbon_road_segment))
mapview::mapview(network_subset)
network_subset_ls = sf::st_cast(network_subset, to = "LINESTRING")
mapview::mapview(network_subset_ls)
network_subset_ls$slopesNASA = slope_raster(network_subset_ls, e = dem_subset) # different error:
Robinlovelace commented 3 years ago

Should be fixed in the latest version:

# install latest version
remotes::install_github("itsleeds/slopes")
library(sf)
library(raster)
library(slopes)

#import datasets
network_original = st_read("https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true")
network = st_zm(network_original, drop = T) # make sure it has no Z values stored
network
u = "https://github.com/U-Shift/Declives-RedeViaria/blob/main/raster/LisboaNASA_clip.tif?raw=true"
download.file(u, "LisboaNASA_clip.tif")
dem = raster("LisboaNASA_clip.tif")
dem

# do they overlap?
raster::plot(dem)
plot(sf::st_geometry(network), add = TRUE)
# why this error?
network$slopeNASA = slope_raster(network, e = dem)
network_ls = sf::st_cast(network, "LINESTRING")
network_ls$slopeNASA = slope_raster(network_ls, e = dem)
plot(network["slopeNASA"], lwd = 5)
Robinlovelace commented 3 years ago

Result (finally after a few changes to the faulty final lines in the code chunk above):

# install latest version
remotes::install_github("itsleeds/slopes")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'slopes' from a github remote, the SHA1 (1877c90f) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(raster)
#> Loading required package: sp
library(slopes)

#import datasets
network_original = st_read("https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true")
#> Reading layer `RedeViariaLisboa_dadosabertos' from data source `https://github.com/U-Shift/Declives-RedeViaria/blob/main/shapefiles/RedeViariaLisboa_dadosabertos.gpkg?raw=true' using driver `GPKG'
#> Simple feature collection with 22098 features and 23 fields
#> geometry type:  MULTILINESTRING
#> dimension:      XYZ
#> bbox:           xmin: -9.229783 ymin: 38.69186 xmax: -9.090931 ymax: 38.79573
#> z_range:        zmin: 0 zmax: 213.8635
#> geographic CRS: WGS 84
network = st_zm(network_original, drop = T) # make sure it has no Z values stored
network
#> Simple feature collection with 22098 features and 23 fields
#> geometry type:  MULTILINESTRING
#> dimension:      XY
#> bbox:           xmin: -9.229783 ymin: 38.69186 xmax: -9.090931 ymax: 38.79573
#> geographic CRS: WGS 84
#> First 10 features:
#>    OBJECTID_1 OBJECTID COD_SIG IDTIPO COD_VIA               DESIGNACAO
#> 1           1    57475  155181      8   67094  Largo de Dona Estefânia
#> 2           2    57430  155182      8   67144  Estrada das Laranjeiras
#> 3           3    57431  155183      8   67144  Estrada das Laranjeiras
#> 4           4    57433  155185      8   67144  Estrada das Laranjeiras
#> 5           5    57434  155186      8   67144  Estrada das Laranjeiras
#> 6           6    57435  155187      8  500000  Rua Interior Particular
#> 7           7    57436  155188      8   67372 Rua Professor Lima Basto
#> 8           8    57437  155189      8   67333         Rua de Campolide
#> 9           9    57438  155190      8  101018                    Via 1
#> 10         10    57439  155191      8   67333         Rua de Campolide
#>    ESTADO_VIA             DTM_ADD             DTM_UPD COMPRIMENT Max_Z Min_Z
#> 1   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      15.58 66.54 66.48
#> 2   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      10.52 76.63 76.13
#> 3   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      31.01 76.13 75.54
#> 4   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      12.90 74.98 74.82
#> 5   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      82.63 75.90 73.37
#> 6   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      17.52 73.37 72.85
#> 7   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      53.49 60.00 58.44
#> 8   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00     103.97 60.56 54.95
#> 9   Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      42.88 57.09 55.65
#> 10  Existente 2012-06-19 01:00:00 2012-06-19 01:00:00      39.27 63.10 62.04
#>    Len_Up Len_Down H_Up H_Down Max_S_Up Max_S_Down Av_S_Up Av_S_Down     slope
#> 1    8.22     7.36 0.06   0.02     0.88       0.36    0.42      0.19 0.3851091
#> 2    0.00    10.54 0.00   0.50     0.00       3.98    0.00      2.72 4.7528517
#> 3    0.00    31.02 0.00   0.59     0.00       3.17    0.00      1.09 1.9026121
#> 4    0.00    12.90 0.00   0.16     0.00       1.55    0.00      0.69 1.2403101
#> 5   23.96    58.87 1.37   2.97     9.11      12.25    3.28      2.89 3.0618419
#> 6   17.53     0.00 0.52   0.00     4.02       0.00    1.70      0.00 2.9680365
#> 7   40.72    12.89 2.01   0.45    14.26       6.09    2.83      2.01 2.9164330
#> 8   35.76    68.77 2.67   5.70    16.23      13.52    4.30      4.77 5.3957872
#> 9   42.91     0.00 1.44   0.00     5.86       0.00    1.92      0.00 3.3582090
#> 10   0.00    39.29 0.00   1.06     0.00       6.53    0.00      1.55 2.6992615
#>                                GlobalID Shape__Length
#> 1  bae5bf0d-ba6f-45ab-b2dd-d2b720c7cbf9      15.57893
#> 2  e071617d-7c01-4a00-9ab8-ee3e9d7486a3      10.52314
#> 3  87a0b0ef-7f92-4426-8d3e-0e4e8cea8d8c      31.01487
#> 4  f95b4446-c6e9-4139-acd7-743c5befecb4      12.89881
#> 5  1c406246-fb27-40d5-80e6-4d31a7134326      82.63346
#> 6  da7a1966-456d-45ed-8347-36690f891871      17.51954
#> 7  ca170299-ed33-48ba-a84d-ee7889f97015      53.48968
#> 8  b2d73e31-149c-47bc-9e65-9a9e00db3915     103.97034
#> 9  c6fcf1c2-15c8-40b1-94b9-9a00498f9b3a      42.88143
#> 10 00776790-0590-4001-a3d8-37e0ed4ecf61      39.26872
#>                              geom
#> 1  MULTILINESTRING ((-9.140667...
#> 2  MULTILINESTRING ((-9.164314...
#> 3  MULTILINESTRING ((-9.164228...
#> 4  MULTILINESTRING ((-9.163915...
#> 5  MULTILINESTRING ((-9.163915...
#> 6  MULTILINESTRING ((-9.163308...
#> 7  MULTILINESTRING ((-9.166655...
#> 8  MULTILINESTRING ((-9.167458...
#> 9  MULTILINESTRING ((-9.167389...
#> 10 MULTILINESTRING ((-9.168217...
u = "https://github.com/U-Shift/Declives-RedeViaria/blob/main/raster/LisboaNASA_clip.tif?raw=true"
download.file(u, "LisboaNASA_clip.tif")
dem = raster("LisboaNASA_clip.tif")
dem
#> class      : RasterLayer 
#> dimensions : 430, 729, 313470  (nrow, ncol, ncell)
#> resolution : 0.0002777778, 0.0002777778  (x, y)
#> extent     : -9.27875, -9.07625, 38.68903, 38.80847  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : /tmp/RtmpASR2Id/reprex2125672f70cf/LisboaNASA_clip.tif 
#> names      : LisboaNASA_clip

# do they overlap?
raster::plot(dem)
plot(sf::st_geometry(network), add = TRUE)

# why this error?
network$slopeNASA = slope_raster(network, e = dem)
#> Error in stop_is_not_linestring(r): Only works with LINESTRINGs. Try converting with sf::st_cast(object, 'LINESTRING') first
network_ls = sf::st_cast(network, "LINESTRING")
#> Warning in st_cast.sf(network, "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant
network_ls$slopeNASA = slope_raster(network_ls, e = dem)
#> [1] TRUE
#> Loading required namespace: geodist
plot(network_ls["slopeNASA"], lwd = 5)

Created on 2020-10-29 by the reprex package (v0.3.0)