Migrating to terra from raster

[scbrown86]

I've started this process, but I'm not 100% sure it's necessary in our case. You can see what I've done here https://github.com/scbrown86/poems/tree/terra_dev. The problem is that now my terra-dev branch is a number of commits behind the poems-main branch and I don't think it can be easily merged (if at all).

I'm just dumping some thoughts below, but happy to discuss.

• The biggest issue with raster is that it was using the rgeos and rgdal packages which are no longer under development and will be removed from CRAN. As of the 18th September 2022 these dependencies have been removed and raster now uses terra for the same functionality. This process started in 2021. See the raster news here.
• We're not really using rasters for anything other than reading in data, population consitency checks and plotting. Everything else (the time consuming part of running the models) under the hood is more or less done with matrices/vectors/arrays.
• terra doesn't readily support parallel processing without wrapping any spatRaster objects. This is easily done, but in my testing we need to keep a wrapped copy of the region$region_raster object that then needs to be unwrapped on each parallel session to make sure the data is consistent with the region. Any spatRaster objects that would be passed on to any parallel functions need to be wrapped. This adds an overhead in wrapping and unwrapping particularly with large regions
• Following that, the poems region object in terra-dev is now initalised with a packed version of the region raster see here
• The package used to do the friction calculations doesn't use terra and still relies on raster for it's calculations. If we swap to terra we would need to convert back to raster before doing those calculations adding additional overhead. At the moment in my terra-dev branch of poems I do it like this no_friction_rast <- raster::raster(terra::unwrap(raster_region$region_raster_packed))
• The friction calculations works fine if calculating frictions sequentially. When the loop is done in parallel it throws an error.
  ## <simpleError in h(simpleError(msg, call)): error in evaluating the argument 'i' in selecting a method for function '[<-': NULL value passed as symbol address>. I'm absolutely sure this is due to the way the conductance raster is stored in the region, but I need to keep working on this when I have time. See here and here
• See Move functions dependecies to terra AgrDataSci/gdistance#3 (comment) for more details on the gdistance issues with swapping to terra
• Without friction the terra-dev version runs without issue which reaffirms to me that the issue is with the conductance raster not being packed.

I'm not overly concerned with raster becoming deprecated particuarly as it's been heavily rewritten to use functionality from terra in place of the old rgeos and rgdal functionality, so maybe it's best to just stick with raster? Particularly given that outside of defining the region, the friction, and plotting, we don't do much else with raster objects?

[japilo]

I have read up on this topic and it seems that we can stay on raster, and probably should, so long as we can lose our dependency on rgeos and rgdal. One obvious way to do this is to set the imported version of raster to be at or after the version where it's changed to no longer depend on these packages. We shall also have to comb through our imported packages and functions to figure out which ones may be causing a dependency on the retiring packages.

[japilo]

I have made a list of imported/suggested packages with problematic dependencies.

• gdistance
• geosphere
• raster v3.4.5 (current version is okay)

So we just need to tick up the version of raster that poems depends on and find replacements for the imported gdistance and geosphere packages.

[scbrown86]

I don't think we need a replacement gdistance as they're slowly adding some additonal methods to handle the conversion from raster to terra. This process started a while ago AgrDataSci/gdistance#3. It means we won't need to convert to terra and then back to raster if the conversion is handled on their end.

The dev version of leastcostpath seems to do something similar but uses terra instead of raster which means we would have to add a function to do the conversion, which possibly (?) brings up the issue of parallel processing with terra objects that I mentioned above. They also don't interally do the conversion if a raster object is passed.

leastcostpath::create_cs(x = as.raster(r), neighbours = 16, dem = NULL, max_slope = NULL)
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘values’ for signature ‘"raster"’

The output of leastcostpath::create_cs() provides a conductance matrix - see the create_cs function and example. I don't know how much effort would be involved in rewriting poems to use this and the output from create_cs. But as it uses terra I think switching to it is more hassle than it's worth.

The geosphere package is written by the author of raster and terra so I would be surprised if he let it become obsolete. Let's assume that it does however - I think the best alternative (and IMO something we really should have done from the start), would be to force only projected coordinates to be used for any spatial model. My understanding is geosphere is only used in poems for calculating dispersal distances on non-projected data? If people want to use poems to make a spatially explicit model, it must use projected coordinates or else the region won't be created and the model won't run. I think the easist way would be through a check in the Region object.

Adding a check to here to make sure the CRS is projected if a region raster is provided: https://github.com/GlobalEcologyLab/poems/blob/b0e18f52b2877593e6f750548f3bdc63d750607f/R/Region.R#L178C1-L178C1
Would force the region to be projected. Something like this.

region_raster = function(value) {
  if (missing(value)) {
    if (is.null(private$.region_raster) && !is.null(private$.coordinates) && self$use_raster) {
      ## if coordinates are provided assume EPSG:4087. Issue warning to use raster object instead
      ## https://epsg.io/4087
      warning("CRS for cooordinates assumed to be EPSG:4087 (https://epsg.io/4087).\nSuggest providing a raster instead of coordinates.",
              immediate. = TRUE, call. = FALSE)
      raster::rasterFromXYZ(cbind(private$.coordinates, cell = 1:nrow(private$.coordinates)),
                            crs = "EPSG:4087")
      
    } else {
      private$.region_raster
    }
  } else {
    if (raster::isLonLat(value)) {
      stop("Region raster must use a projected CRS", call. = FALSE)
    }
    if (!is.null(value) && !("RasterLayer" %in% class(value))) {
      stop("Region raster should be a raster::RasterLayer (or inherited class) object", call. = FALSE)
    }
    if (!is.null(value) && !is.null(private$.coordinates)) {
      stop("Region is already associated with a set of coordinates", call. = FALSE)
    }
    if (raster::fromDisk(value)) { # move to memory
      raster::values(value) <- raster::values(value)
    }
    private$.region_raster <- value
  }
}

We would need an additional check if coordinates instead of a raster is passed to the region. Something to stop if the x coordinates are within -365 and 365 and the y coordinates are within -90.1 and 90.1.

[japilo]

So in summary, it seems the problem will resolve on its own if we just wait a bit for geosphere and gdistance to get off of rgeos and rgdal.

[japilo]

The new plan is to keep raster for now and shift functions away from geosphere and gdistance, which are still stuck on using sp.

[scbrown86]

Like I mentioned in my earlier comment, the geosphere dependency can be resolved by forcing people/groups to use projected coordinates. It is only used for calculating distances when lat/long rasters or coordinates are given. It's a really simple fix, and in my opinion, a conditon which probably should've been enforced from the beginning of package development.

There is also no need to migrate away from gdistance as sp isn't being removed from CRAN or even retired. Internally, sp has been migrated to use sf instead of rgdal, rgeos, and maptools. Those three packages are being retired and removed very soon, any packages which don't require a sp version after those functions have been transitioned to use sf will fail. gdistance only uses a single sp function which will be affected by the transition - sp::CRS. See here - https://r-spatial.org/r/2023/05/15/evolution4.html. As long as the latest version of sf is used, sp will use the sf package to extract the CRS and it won't be an issue. The easist fix is to force a version dependency for sp and sf in poems DESCRIPTION file - probably under suggests as we don't actually use the packages in any functions, but it will force any downstream packages to use the latest versions.

From a quick look, gdistance only uses sp::CRS when coordinates are given or for showing the least cost path between two points. The transition layers, barriers, cost distance etc. are all created from raster inputs. As poems passes a raster through to all the gdistance functions, poems is not affected by any of the changes.

[damienfordham]

Sounds like the right way forward Stu. From: Stuart Brown ***@***.***> Date: Sunday, 24 September 2023 at 10:19 am To: GlobalEcologyLab/poems ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [GlobalEcologyLab/poems] Migrating to terra from raster (Issue #4) CAUTION: External email. Only click on links or open attachments from trusted senders.

…

________________________________ Like I mentioned in my earlier comment, the geosphere dependency can be resolved by forcing people/groups to use projected coordinates. It is only used for calculating distances when lat/long rasters or coordinates are given. It's a really simple fix, and in my opinion, a conditon which probably should've been enforced from the beginning of package development. There is also no need to migrate away from gdistance as sp isn't being removed from CRAN or even retired. Internally, sp has been migrated to use sf instead of rgdal, rgeos, and maptools. Those three packages are being retired and removed very soon, any packages which don't require a sp version after those functions have been transitioned to use sf will fail. gdistance only uses a single sp function which will be affected by the transition - sp::CRS. See here - https://r-spatial.org/r/2023/05/15/evolution4.html<https://r-spatial.org/r/2023/05/15/evolution4.html>. As long as the latest version of sf is used, sp will use the sf package to extract the CRS and it won't be an issue. The easist fix is to force a version dependency for sp and sf in poems DESCRIPTION file - probably under suggests as we don't actually use the packages in any functions, but it will force any downstream packages to use the latest versions. From a quick look, gdistance only uses sp::CRS when coordinates are given<https://github.com/AgrDataSci/gdistance/blob/cef36d7c27ec2ff5ab0a3e072afe21328c06de29/R/internal-functions.R#L27> or for showing the least cost path between two points<https://github.com/AgrDataSci/gdistance/blob/cef36d7c27ec2ff5ab0a3e072afe21328c06de29/R/shortestPath.R#L140>. The transition layers, barriers, cost distance etc. are all created from raster inputs. As poems passes a raster<https://github.com/GlobalEcologyLab/poems/blob/b1964bd16979f9aafaea73fda1baf13fd918f227/R/DispersalFriction.R#L118C9-L118C9> through to all the gdistance functions, poems is not affected by any of the changes. — Reply to this email directly, view it on GitHub<#4 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AEDXBYCB6CBVKKZ54K5S75LX357QTANCNFSM6AAAAAA3MJXKCI>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>