Compare with ArcMap method

[Robinlovelace]

In terms of:

• [x] Similarity with 'Avg_Slope'
• [ ] Other measures of slope
• [x] Speed of execution

[Robinlovelace]

Timing on 3410 routes:

# routes = readRDS("~/wip/pctLisbon-data/routes_cyclestreet_segments_lisbon.Rds")
r = sf::read_sf("~/wip/pctLisbon-data/route_segments_gradient2_projected.gpkg")
e = raster::raster("~/wip/pctLisbon-data/Shapefiles/DEM_ist_10m/r1/")
nrow(routes)
#> Error in nrow(routes): object 'routes' not found
system.time({
  s = slopes::slope_raster(r, e)
})
#>    user  system elapsed 
#>   1.142   0.020   1.169

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

[temospena]

In ArcGIS 10.1 SP1, with 3D extension

3 seconds (4641 routes) (I projected in EPSG:3763 previously, to be compatible with the DEM)

[Robinlovelace]

Could you try running the R code on the same computer so it's a fair comparison?

[temospena]

good point. yes...

On Mon, May 4, 2020 at 7:56 PM Robin notifications@github.com wrote:

Could you try running the R code on the same computer so it's a fair comparison?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/ITSLeeds/slopes/issues/10#issuecomment-623643100, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJKLUXSP4XF3SONI46QFFVTRP4FU7ANCNFSM4MY6PRWQ .

[temospena]

After all the other one had 4641 routes!

r = sf::read_sf("D:/GIS/Lovelace/route_segments_gradient2_projected/lisbon_road_segments_TM06.shp")
> e = raster::raster("D:/GIS/Lovelace/pctLisbon-data/Shapefiles/DEM_ist_10m/r1/w001001.adf")
> nrow(r)
[1] 4641
> system.time({
+   s = slopes::slope_raster(r, e)
+ })
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 09s
   user  system elapsed 
   7.94    2.37   10.34 
> 

so, ArcGIS (4641 routes) --- 3 sec slopes(4641 routes) ---10 sec

I can try again with the 3410 routes, but let's try with a bigger dataset

[temospena]

route_segments_gradient2_projected.gpkg")

for the same dataset with 3410 routes 2 sec.

[Robinlovelace]

OK, good motivation for speed-ups with terra package :rocket:

[Robinlovelace]

Here are the addition columns from ArcMap that would be worth exploring to see if we can reproduce the results. Good basis for description in the vignette also:

library(slopes)
summary(lisbon_road_segments[4:11])
#>    Shape_Leng          Z_Min            Z_Max            Z_Mean      
#>  Min.   :   1.81   Min.   : 1.225   Min.   : 3.688   Min.   : 2.508  
#>  1st Qu.:  39.52   1st Qu.: 5.012   1st Qu.: 9.196   1st Qu.: 6.475  
#>  Median :  95.24   Median :14.372   Median :16.032   Median :15.000  
#>  Mean   : 237.89   Mean   :26.975   Mean   :31.405   Mean   :29.207  
#>  3rd Qu.: 198.85   3rd Qu.:46.242   3rd Qu.:52.697   3rd Qu.:48.598  
#>  Max.   :3749.71   Max.   :94.280   Max.   :95.019   Max.   :94.704  
#>     SLength           Min_Slope          Max_Slope        Avg_Slope      
#>  Min.   :   1.834   Min.   : 0.00000   Min.   : 0.000   Min.   : 0.0000  
#>  1st Qu.:  39.528   1st Qu.: 0.02576   1st Qu.: 2.188   1st Qu.: 0.9728  
#>  Median :  90.722   Median : 0.35968   Median : 5.769   Median : 2.7073  
#>  Mean   : 236.913   Mean   : 1.69648   Mean   : 8.869   Mean   : 4.2789  
#>  3rd Qu.: 191.581   3rd Qu.: 1.71447   3rd Qu.:13.720   3rd Qu.: 6.4613  
#>  Max.   :3751.715   Max.   :21.48414   Max.   :46.021   Max.   :21.4841

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

[temospena]

Timing with about 25k segments

ArcGIS (24927 routes) --- 11 sec slopes(24927 routes) --- 98 sec

But 100 seconds is good anyway to process a whole city.

RedeViaria = read_sf("~/Declives_Lisboa/Declives_RedeViariaLisboa_TM06_RF.shp")
st_crs(RedeViaria) = 3763 #TM06, same as r1
nrow(RedeViaria) #24927
system.time({
  sRedeViaria25k = slopes::slope_raster(RedeViaria, e)
})
#user  system elapsed 
#71.88   26.04   98.72

[Robinlovelace]

Thanks Rosa, just to check, is that using the latest version of the package installed with:

devtools::install_github("itsleeds/slopes")

?

Imagine so. I think switching to terra will speed it up, we can keep this issue open and document any performance things here.

[Robinlovelace]

Heads-up, I get this based on a slightly different dataset with ~4600 route segments:

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

# aim: test performance
r = read_sf("~/wip/pctLisbon-data/Shapefiles/lisbon_segments_red_Zelevation.shp")
e = raster::raster("~/wip/pctLisbon-data/Shapefiles/DEM_ist_10m/r1")

nrow(r)
#> [1] 4641
# st_crs(r) = 3763 #TM06, same as r1
r = sf::st_transform(r, sf::st_crs(lisbon_road_segment))
nrow(r) # 4641
#> [1] 4641
system.time({
  s = slopes::slope_raster(r, e)
})
#>    user  system elapsed 
#>   1.876   0.104   1.999

summary(s)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#> 0.00000 0.01235 0.02608 0.03568 0.04776 0.27673       2

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

[temospena]

After updating slopes, here it goes for that dataset

##Again with 4600 routes
library(slopes)
library(sf)
r4600 = read_sf("~/Shapefiles/lisbon_segments_red_Zelevation.shp")
nrow(r4600) #4641
e = raster::raster("~/Shapefiles/DEM_ist_10m/r1/w001001.adf")

r4600 = sf::st_transform(r4600, sf::st_crs(lisbon_road_segments))
nrow(r4600) # 4641
system.time({
  s4600 = slopes::slope_raster(r4600, e)
})
#user  system elapsed 
#2.78    0.19    2.98 
summary(s4600)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#0.00000 0.01235 0.02608 0.03568 0.04776 0.27673       2 

[Robinlovelace]

And in ArcMap? My preliminary interpretation is that our package is comparable speed and accuracy as an established proprietary solution, and we have not yet optimised it using the terra package as a backend.

[temospena]

sorry, in ArcMap --- 3sec also!

[Robinlovelace]

Muy interesante!

[temospena]

And again, after updating slopes For the 25k dataset of segments (not entire routes) It reduced from 98sec to 67 sec

RedeViaria = read_sf("~/Declives_RedeViariaLisboa_TM06_RF.shp")
r25k = sf::st_transform(RedeViaria, sf::st_crs(lisbon_road_segments))
st_crs(r25k) #TM06, 3763
nrow(r25k) #24927
system.time({
  s25k = slopes::slope_raster(r25k, e)
})
#user  system elapsed 
# 55.57   11.48   67.08 
summary(s25k)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  0.000   0.000   0.000   0.002   0.000   0.015   24917 

But having an issue with NAs...

(it happened the same with sRedeViaria25k)

[temospena]

Good news!

I discarded some road segments from the 25k dataset (the ones that were pretty close to the DEM border, i.e, the river), and kept about 23k. Updated the package.

Run in Arc and in slopes pkg ArcGIS (23k routes) --- 11 sec slopes(23k routes) ---70 sec

Here are the performance results

RedeViaria23 = read_sf("lisbon_road_segments_TM06_23k.shp")
r23k = sf::st_transform(RedeViaria23, sf::st_crs(lisbon_road_segments))
st_crs(r23k) #TM06, 3763
e = raster::raster("~/Shapefiles/DEM_ist_10m/r1/w001001.adf")
nrow(r23k) #23058
system.time({
  s23k = slopes::slope_raster(r23k, e)
})
#  user  system elapsed 
#50.32   19.70   70.17 
summary(s23k)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.00000 0.01443 0.03163 0.04429 0.05805 0.79608 

No NAs this time!

And here is the correlation value

r23k$slope = s23k
summary(r23k$slope*100) #slopes pkg
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.000   1.443   3.163   4.429   5.805  79.608 
summary(r23k$Avg_Slope) #ArcGIS
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.000   1.444   3.165   4.440   5.807  79.608 

cor(r23k$slope*100, r23k$Avg_Slope)
#0.999914

:fireworks:

[Robinlovelace]

That is an impressive result @temospena! I think we can see this issue as 'closed' when these results are reported in the 'validation' vignette.

[Robinlovelace]

P.s. do you mind if I tweet about this? Can tag you if you have an account. Want to tell the world!

[temospena]

sure, no problem, it's @rosamfelix

[Robinlovelace]

Shape_Leng Z_Min Z_Max Z_Mean Are the columns from ArcMap