landscitech / roads Goto Github PK

What should happen with these? I think this typically means that an area should not be used for roads, in which case the NA should be replaced with a very high cost. There is currently a error when there is an NA cost in a road path. This happens when a road that is burned in goes through an NA pixel because it will stay NA after multiplying the road and cost rasters. But should probably warn if any NAs detected in cost raster. But how time consuming is it to check all not NA?

Files containing functions derived from Kyle's code should be prefaced with copyright/licensing notice. See getGraph for example.
All functions derived from Kyle's code were named roadCLUS., and contained in helperFnsClus.R file. Now called getGraph, lcpList, mstList, shortestPaths, getClosestRoad, buildSnapRoads.

Other R files should also include header, with her majesty 2021 as copyright holder.

Got this error on Fedora Linux, R-devel, clang, gfortran platform

4432#> * installing *source* package ‘s2’ ...

4433#> ** package ‘s2’ successfully unpacked and MD5 sums checked

4434#> ** using staged installation

4435#> Testing compiler using PKG_CFLAGS=

4436#> --------------------------- [ANTICONF] --------------------------------

4437#> Configuration failed because openssl was not found. Try installing:

4438#> * deb: libssl-dev (Debian, Ubuntu, etc)

4439#> * rpm: openssl-devel (Fedora, CentOS, RHEL)

4440#> * csw: libssl_dev (Solaris)

4441#> * brew: [email protected] (Mac OSX)

4442#> If openssl is already installed, check that 'pkg-config' is in your

4443#> PATH and PKG_CONFIG_PATH contains a openssl.pc file. If pkg-config

4444#> is unavailable you can set INCLUDE_DIR and LIB_DIR manually via:

4445#> R CMD INSTALL --configure-vars='INCLUDE_DIR=... LIB_DIR=...'

4446#> -------------------------- [ERROR MESSAGE] ---------------------------

4447#> tools/version.c:1:10: fatal error: 'openssl/opensslv.h' file not found

4448#> #include

4449#> ^~~~~~~~~~~~~~~~~~~~

4450#> 1 error generated.

4451#> --------------------------------------------------------------------

4452#> ERROR: configuration failed for package ‘s2’

4453#> * removing ‘/home/docker/R/s2’

Val thinks it will be fine though

Hi Sarah,

While it is on my mind, I just wrote a method for getting distance from roads (or any sources) for David that is a lot faster (seconds not days) than the method he has been using. See details below. I think we should add this to the roads package, as it will likely come up again. However, it uses our parallel focal method/package that Val is polishing now. We need to release that package before adding as an optional dependency to roads.

Ok, time for me to stop fiddling with fun methods and go to bed.

Onward,

Josie
From: Hughes,Josie (ECCC)
Sent: August 5, 2021 10:32 PM
To: Lapins,David (ECCC) [email protected]
Subject: RE: Revelstoke Binary Rasters for Distance Metric

Ok that was entertaining.

I have added a getDistanceFromSource.R to scripts folder in the RoadAnalysis github repository.

The function uses pfocal, which is a fast parallel moving window function we have been developing for another project.
In order to install the pfocal package, you will need to make sure you have a current version of R installed. If you run into problems installing (and you may) let us know. Now is a good time for us to find and fix problems with that package.
There is a tradeoff between speed and accuracy in this method, controlled by the parameters kwidth and dissag.
kwidth is the radius of the moving window in number of cells. For very small values (e.g. 1,2) we get grid artefacts in the results. Larger moving windows are more circular which reduces the grid artefacts. If dissag is false, the resulting distance map will have values of 0, 1kwidth, 2kwidth, etc. If dissag is true the algorithm disaggregates the raster to a finer resolution, so the resulting distance map will have values of 0,1,2, etc, but there is considerable computational cost associated with the higher accuracy.

On my laptop with 4 cores the big map calculation at the end of the script takes 36 seconds. I will be interested to hear how much slower it is on your machine.

You should probably do some comparisons to your other distance results for the small map. This looks roughly right but I haven’t done any formal testing.

Let me know how it goes.

best,
Josie

In both getLandingsFromTarget and projectRoads I get some weirdness when st_geometry_type(landings, by_geometry = FALSE) returns "GEOMETRY" Should be an error that says geotype can't be GEOMETRY and try st_cast

TODO: Make this and a list of sf objects, (or maybe a year column?) work

 # ### using: scenario 3 / RasterStack landings / minimum spanning tree ("mst")
 # # demo scenario 3
scen <- demoScen[[3]]

 # # landing sets 1 to 4 of scenario 3, as RasterStack:
land.rstack <- scen$landings.stack[[1:4]]

prRes <- projectRoads(land.rstack, scen$cost.rast, scen$road.rast ,"mst",
                         plotRoads = doPlots, mainTitle = "Scen 3: RasterStack-MST")

When checking "roads" with the dev version of "raster" (https://github.com/rspatial/raster) I now get

  ══ Failed tests ════════════════════════════════════════════════════════════════
  ── Failure ('test-errors-messages.R:8:3'): different CRS is error ──────────────
  `projectRoads(CLUSexample$landings, cost, rds)` did not throw an error.

This is because of changes in "raster" related to the changes in the "sp family" packages. sp::identicalCRS does not work well with Raster* objects any more.

You now need to use CRS(Raster*) or you could use terra::same.crs(x, y). (The best approach might be to drop "raster" and only use "terra").

Can you please fix this so that your package won't break when I submit "raster" to CRAN?

In the roads manuscript it says: "Each STAP (single target access problems) was then solved sequentially, with the newly projected roadway updating the existing road network and subsequently the costs of road construction." This is not actually happening in the package when "lcp" method is used. The cost is updated after the shortest path to each landing is calculated so if all the landings are passed in one call to projectRoads then the updated costs are not used. We can use updated costs if the object produced in the first call to projectRoads is supplied to the function as the sim argument. This has been the case since Kyles example code here. Not sure if this is something we would want to change in the package or if we just need to make it more clear in the paper that we used the least cost path version where each landing is considered independently of others.

Here is an example of how we do it now vs how it would be done with cost updating after each landing is built.

using the master branch of roads

# Check on the LCP method make sure that it is actually dynamically updating
# cost after each individual road is built

library(roads)
library(dplyr)
library(raster)

rastColours2 <- c('grey50', viridis::viridis(30))

## scenario 
scen <- demoScen[[1]]

## landing sets 2 
land.pnts2 <- scen$landings.points[scen$landings.points$set == 2,]

oneTime_sim <- projectRoads(land.pnts2,
                                   scen$cost.rast,
                                   scen$cost.rast==0,
                                   roadMethod='lcp')$roads


## plot
oneTime_sim <- subs(oneTime_sim, data.frame(from = c(0, 1), to = c(NA,1)))
plot(scen$cost.rast, col = rastColours2, 
     main = "Cost not updated after each road built",
     legend = FALSE)
plot(oneTime_sim, add = TRUE, col = "grey50", legend = FALSE)
plot(land.pnts2, add = TRUE,
     pch = 21, cex = 1.5, bg = 'white')
text(land.pnts2@coords, labels = land.pnts2$ID, cex = 0.6, adj = c(0.5, 0.3),
     xpd = TRUE)

# doesn't look like cost is dynamically updated

# Try iterating over individual landings b/c that should definitely update it

iterLands_sim <- list(projectRoads(land.pnts2[land.pnts2$ID==1,],
                                   scen$cost.rast,
                                   scen$cost.rast==0,
                                   roadMethod='lcp', roadsOut = "raster")) 
for (i in 2:max(land.pnts2$ID)){
  iterLands_sim <- c(iterLands_sim,
                   list(projectRoads(sim = iterLands_sim[[i-1]], 
                                landings = land.pnts2[land.pnts2$ID==i,], 
                                roadsOut = "raster")))
}

## plot
plot(scen$cost.rast, col = rastColours2, 
     main = "Update cost for each landing",
     legend = FALSE)
iterLands_sim[[7]]$roads <- subs(iterLands_sim[[7]]$roads, 
                    data.frame(from = c(0, 1), to = c(NA,1)))
plot(iterLands_sim[[7]]$roads, add = TRUE, col = "grey50", legend = FALSE, lwd = 2)
plot(land.pnts2, add = TRUE,
     pch = 21, cex = 1.5, bg = 'white')
text(land.pnts2@coords, labels = land.pnts2$ID, cex = 0.6, adj = c(0.5, 0.3),
     xpd = TRUE)

Old version of getLandingsFromTarget

rast <- demoScen[[1]]$landings.poly %>% raster::rasterize(demoScen[[1]]$cost.rast)
  # make sure that a single celled havest block will work with clumps
  rast[10,10] <- 6
  outOld <- getLandingsFromTarget(rast > 0, numLandings = 614*0.1)

New version

  rast <- demoScen[[1]]$landings.poly %>% raster::rasterize(demoScen[[1]]$cost.rast)
  # make sure that a single celled havest block will work with clumps
  rast[10,10] <- 6
  
  out <- getLandingsFromHarvest(rast, 0.25)

The new version is quite different but I think is an improvement since it uses a regular sampling grid to get points with in the polygon rather than random sampling so the resulting points are spread out in the polygon rather than some of them ending up side by side which isn't very realistic. It is also easier to set numLandings since it is the density of landings

It took 3.5 mins to get 90777 points in 10338 polygons takes longer the higher the density of points

This isn't so much an issue as a request for comment. If @josie-hughes you think it is ok I will rename it to getLandingsFromTarget and remove the old version. Old code will have to be changed to reflect the new meaning of numLandings which is 1/[distance between points on the grid] (we could make it simply the distance between points in the grid but I was thinking of it as density so smaller equals further apart)

pfocal should be in Suggests but that doesn't work when it is not on CRAN so it should be added later. The current plan allows pkgdown to work but fails to pass the check.

We only use one fun SpaDES.tools::adj and there is a terra version that might work. So compare and then change

Right now we do not change the cost in the raster when new roads are built. This is helpful when plotting the result because you can see the cost that was used to project the roads. However it slows things down when the output roads are in raster format and when a sim list is passed to projectRoads in order to use the existing graph (which is updated as roads are built). If the cost surface was updated for new roads then we could output the cost == 0 as the raster roads layer and we would not have to burn in the roads from the previous simulation when a sim list is used as input.
@josie-hughes do you think that would be a good change?

First release:

• usethis::use_cran_comments()
• Proof read Title: and Description:
• Check that all exported functions have @returns and @examples
• Check that Authors@R: includes a copyright holder (role 'cph')
• Check licensing of included files
• Review https://github.com/DavisVaughan/extrachecks

Prepare for release:

• devtools::build_readme()
• urlchecker::url_check()
• #12
• devtools::check_win_devel()
• #13
• Draft blog post

Submit to CRAN:

• usethis::use_version('minor')
• devtools::submit_cran()
• Approve email

Wait for CRAN...

• Accepted 🎉
• usethis::use_github_release()
• usethis::use_news_md()
• usethis::use_dev_version()
• Update install instructions in README
• Finish blog post
• Tweet
• Add link to blog post in pkgdown news menu

see https://github.com/LandSciTech/RoadAnalysis/blob/master/scripts/getDistanceFromSource.R

I have noticed that the results for the CLUS example data are different for Kyle's version of the package compared to our version based on comparing the results in the vignette to kylesCLUSExample.RMD. (This is different from the version here that uses different input data). This difference existed in the version that came from Pat Kirby so I don't think it is a result of the changes I made.

LCP

My version

Pat's version

Kyle's version

MST

My version

Pat's version

Kyle's version

This works for me with the data Josie sent

devtools::load_all(".")
outObj <- readRDS("../../../Downloads/roadInputs.RDS")
roadsCH <- raster::raster("../../../Downloads/churchillRoadRaster.tif")
raster::crs(roadsCH) <- raster::crs(outObj$roads)
raster::extent(roadsCH) <- raster::extent(outObj$roads)
raster::origin(roadsCH) <- raster::origin(outObj$roads)
raster::res(roadsCH) <- raster::res(outObj$roads)
outObj$roads <- roadsCH
sim_mst = projectRoadsNew(outObj$landings,
                      outObj$cost,
                      outObj$roads,
                      roadMethod="mst",
                      roadsOut = "raster")

the roadsOut part shouldn't change anything because it is the default

I get results but they don't look very good with some points not getting connected to the road network

I have noticed though that the cost is 0 in some areas that don't have roads yet which is probably messing things up.
plot(outObj$roads)

plot(outObj$cost == 0)

cost=cMap;cost[cost>0]=0;cost[cMap==0]=100

rds = raster::rasterToPoints(cMap, fun = function(x){x > 0},
                             spatial = TRUE) 
rds=spsample(rds,n=1,type="random")

nw = roads::projectRoadsNew(landings = cMap,
                                   cost = cost,
                                   roads = rds,roadsOut="sf")

I removed the simpleCost function since we weren't using it any where and it is probably too simple to export as part of the package. Should we add a function to facilitate calculating costs though? Or is this a whole other problem best left to the users?
@josie-hughes any thoughts?

Having trouble with the manual says it is having latex issues. Not sure if this will be a problem or just that my system doesn't work but CRANs will?

Start a new branch to work on using terra instead of raster. Continue supporting raster inputs though. Look into whether terra can do the things we currently do with stars to see if we can remove that dependency. I will document any challenges that come up here as a reference for repeating the process in other packages

The terra website has a table showing new function names and raster equivalents and highlights differences in behaviour https://rspatial.github.io/terra/reference/terra-package.html#xxv-new-method-names-1

See sf issue 1730

I changed the vignette to have a slanted line to get around this but we should change it back after.

First release:

• usethis::use_cran_comments()
• Update (aspirational) install instructions in README
• Proofread Title: and Description:
• Check that all exported functions have @return and @examples
• Check that Authors@R: includes a copyright holder (role 'cph')
• Check licensing of included files
• Review https://github.com/DavisVaughan/extrachecks

Prepare for release:

• git pull
• Check if any deprecation processes should be advanced, as described in Gradual deprecation
• devtools::build_readme()
• urlchecker::url_check()
• devtools::check(remote = TRUE, manual = TRUE)
• devtools::check_win_devel()
• rhub::check_for_cran()
• git push
• Draft blog post

Submit to CRAN:

• usethis::use_version('minor')
• devtools::submit_cran()
• Approve email

Wait for CRAN...

• Accepted 🎉
• git push
• usethis::use_github_release()
• usethis::use_dev_version()
• git push
• Finish blog post
• Tweet
• Add link to blog post in pkgdown news menu

This is top left corner of churchillExample.R results. Does a long straight line crossing existing road to connect the two clumps.
this is probably related to the conn and keep part of pathsToLine