jeffreyevans/spatialeco

R package for spatial analysis and modelling of ecological systems

spatialeco's Issues

NA handling when some points are not in polygons

I came across the same problem as mentionned in this post. Looking at the new version of the function, it appears that NA are deleted in statement z <- pts[!is.na(sp::over(pts, sp::geometry(polys))), ] but aren't in the next statement z@data <- data.frame(sp::over(pts, polys)). As a consequence data.frame(sp::over(pts, polys)) and z@data don't have the same number of rows, triggering an error.
In my case NA values comes from the fact that some points do not fall in any of the polygons tested. This a not a bug since I'm testing a whole bunch of points to find which ones fall in the polygons.

I changed the two statement adding a NA removal line that seems to do the trick, but maybe you did that for a reason ?

z <- x[!is.na(sp::over(x, sp::geometry(y))), ]
zover <-  data.frame(sp::over(x, y))
z@data <- data.frame(zover[!is.na(zover[[poly.id]]), ])

cant find package sf

Hi

I am trying to install spatialEco in R running inside a Linux docker container. The installation says that it was successful but spatialEco complains that it cannot find sf. I tried to install sf specifically as well. sf installs fine but I can't still use spatialEco. Here is the error I get:
<packageNotFoundError in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): there is no package called ‘sf’>

Not sure if that is a bug. A sample Dockerfile showing that the install was successful is:

FROM rstudio/plumber
RUN R -e "install.packages('sf')"
RUN R -e "install.packages('spatialEco')"

EXPOSE 8000

On windows, it works fine.

Github repository not mentioned in package description

Dear Jeff,

you might want to add a link to your github repo to the DESCRIPTION file so possible issues can be addressed here.

This line in DESCRIPTION will add a link in the CRAN man page:
BugReports: https://github.com/jeffreyevans/spatialEco/issues

No raster.modifed.ttest function

First of all, thank you for your work. I've seen your name a ton on SE and appreciate your input.

I've used the spatialEco package before. I'm surprised this isn't working. It just says: Error in raster.modifed.ttest(G1, G2) : could not find function "raster.modifed.ttest"

I've tried on two machines, windows and linux. Maybe there is a dependency I'm unaware of or an update I'm missing? Or maybe I'm just making a silly error somewhere. Here is my output from your example code:

 library(gstat)                                         
Warning message:
package ‘gstat’ was built under R version 4.1.2 
> library(sp)
> library(raster)
> library(spatialEco)                                            

Attaching package: ‘spatialEco’

The following object is masked from ‘package:raster’:

    shift

> data(meuse)                                            
> View(meuse)
> data(meuse.grid)                                       
> View(meuse.grid)
> coordinates(meuse) <- ~x + y                           
> coordinates(meuse.grid) <- ~x + y 
> # GRID-1 log(copper):                                              
> v1 <- variogram(log(copper) ~ 1, meuse)                  
> x1 <- fit.variogram(v1, vgm(1, "Sph", 800, 1))           
> G1 <- krige(zinc ~ 1, meuse, meuse.grid, x1, nmax = 30)
[using ordinary kriging]
> gridded(G1) <- TRUE                                      
> G1@data = as.data.frame(G1@data[,-2])
> # GRID-2 log(elev):                                              
> v2 <- variogram(log(elev) ~ 1, meuse)                  
> x2 <- fit.variogram(v2, vgm(.1, "Sph", 1000, .6))        
> G2 <- krige(elev ~ 1, meuse, meuse.grid, x2, nmax = 30)
[using ordinary kriging]
> gridded(G2) <- TRUE    
> G2@data <- as.data.frame(G2@data[,-2])
> G2@data[,1] <- G2@data[,1]
> corr <- raster.modifed.ttest(G1, G2)      
Error in raster.modifed.ttest(G1, G2) : 
  could not find function "raster.modifed.ttest"
> 

Here is the package version installed:

packageVersion("spatialEco")
[1] ‘1.3.7’

Here is session info for the linux system I used in checking this issue if this helps:

sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Thanks again for your presence in the FOSS community and thank you in advance for your insights into this issue. If you need any further info I'm happy to provide.

could not find function "sg.smooth"

I get this error: Error in sg.smooth(t(... : could not find function "sg.smooth" when I run my code using Windows computer. The "spatialEco" package is installed.. Checked the version and it says 1.1-0. I originally wrote my code in my Mac and I didn't experience any glitch. I noticed the package version in my Mac and it says 1.2-0. Could this be causing the problem? How come they have different versions when I just did the same process of install.packages()?

Thanks!

Weighted KDE error

If I'm understanding correctly, the option to pass values to the argument y in sp.kde would give a density plot where each point would have a predicted value based on the values y. If so this is exactly what I need and would be so grateful to make it work. I am hoping to generate a kde plot for air pollution levels based on sampled locations.

I have no problems running unweighted kde passing values to x, nr, and nc. However, as soon as I include values for y, I get an error:

Error in matrix(rep(w, n[1]), nrow = n[1], ncol = nx, byrow = TRUE) * :
non-numeric argument to binary operator

I get this same error even if I remove nr and nc parameters and let the function pick default values.

Help much appreciated.

Problem with bandwidth calculation in sp.kde

Dear Jeff,

in the CRAN version of spatialECO is a bug in the bandwidth calculation for the sp.kde function: the estimation function doesn't make use of the coordinate pairs.

This bug is already fixed in commit 1603c1e.

spatialEco time-series data

Hi Jeffrey,

Currently, I am working on time-series data. I have created raster stack with LST as dependent variable (each band corresponds to acquisition date) and raster stack with 2 explanatory variables (NDVI and DEM saved as raster stack in list) for every date as well. Now, I am trying to use loop to downscale my time-series images.

However, after iteration of 75 layers (I have 225 dates) I get the error: ''Error in seq_len(p) : argument must be coercible to non-negative integer''. The problem is strange because when I try to downscale the image (as 76th layer) separately I don't have problem with downscaling.

I would be grateful for help.

Best,
Paulina

breeding.density - "implicit list embedding of S4 objects is deprecated"

Hi @jeffreyevans,

Looks like changes to R as of version 3.3.0 introduced an annoying warning to breeding.density. Minor, but looks to be a one-liner fix.

Best,
Jason

Reproduce issue
The issue can be reproduced using the help doc example for breeding.density (I'm using R version 3.5.2 and platform i386-w64-mingw32/i386 (32-bit)). The returned object appears correct, but >= 50 warnings are raised:

In `[<-`(`*tmp*`, i, value = sp::Polygons(list(polys[[i]]),  ... :
  implicit list embedding of S4 objects is deprecated

Cause
The warning appears to be caused by a change to how R handles S4 objects beginning in R version 3.3.0. Link to more info.

Solution
A bit of testing suggests the warning can be removed by modifying line 83 of the function (as of spatialEco 1.1.1) to use double brackets:

Single brackets, as in line 83 of breeding.density in spatialEco version 1.1.1
pl[i] <- sp::Polygons(list(polys[[i]]), i)

Double brackets, modified to comply with R>3.3.0 S4 standards, eliminates warning
pl[[i]] <- sp::Polygons(list(polys[[i]]), i)

Calculation of global spatial cross-correlation index (GSCI)

Hey!

I was wondering about the calculation of GSCI, in particular line 157 of the corresponding R code of function crossCorrelation (see below).
Where does the second part stem from [ / ( length(x) - 1 ) ]?
As far as I can tell, this term is not a part of the formulas shown in Chen (2015). I also tested the function with the original data given by Chen (2015). The GSCI given by the function was 0.005594442, the GSCI given in Chen (2015) is 0.1566. This is in fact 0.005594442 * (length(x) - 1), which is 0.005594442 * (29 - 1).

global.i <- as.numeric(x%%w%%y) / (length(x) - 1)

Install error

Hi,

I am trying to install spatialEco package but I receive the following error message. I am working in Ubuntu with all the updates up-to-date.

• installing source package ‘sf’ ...
  ** package ‘sf’ successfully unpacked and MD5 sums checked
  ** using staged installation
  configure: CC: gcc -std=gnu99
  configure: CXX: g++ -std=gnu++11
  checking for gdal-config... no
  no
  configure: error: gdal-config not found or not executable.
  ERROR: configuration failed for package ‘sf’
• removing ‘/home/jose/R/x86_64-pc-linux-gnu-library/4.0/sf’
  Warning in install.packages :
  installation of package ‘sf’ had non-zero exit status
• installing source package ‘rgeos’ ...
  ** package ‘rgeos’ successfully unpacked and MD5 sums checked
  ** using staged installation
  configure: CC: gcc -std=gnu99
  configure: CXX: g++ -std=gnu++11
  configure: rgeos: 0.5-5
  checking for /usr/bin/svnversion... no
  configure: svn revision: 640
  checking for geos-config... no
  no
  configure: error: geos-config not found or not executable.
  ERROR: configuration failed for package ‘rgeos’
• removing ‘/home/jose/R/x86_64-pc-linux-gnu-library/4.0/rgeos’
  Warning in install.packages :
  installation of package ‘rgeos’ had non-zero exit status
  ERROR: dependency ‘sf’ is not available for package ‘spdep’
• removing ‘/home/jose/R/x86_64-pc-linux-gnu-library/4.0/spdep’
  Warning in install.packages :
  installation of package ‘spdep’ had non-zero exit status
  ERROR: dependencies ‘sf’, ‘spdep’, ‘rgeos’ are not available for package ‘spatialEco’
• removing ‘/home/jose/R/x86_64-pc-linux-gnu-library/4.0/spatialEco’
  Warning in install.packages :
  installation of package ‘spatialEco’ had non-zero exit status

The downloaded source packages are in
‘/tmp/RtmpRUuLIk/downloaded_packages’

Thank you beforehand,
Jose Miguel

crossCorrelation

It seems global.i and global.p are wrong. They used to be defined as (https://github.com/jeffreyevans/spatialEco/blob/42641f1429160627e4459a2bfad24d114a61e5bc/R/crossCorrelation.R):

global.i <- as.numeric(x%*%w%*%y) / (length(x) - 1)
( global.p <- sum( abs(isim) > abs( global.i ) ) / length(isim) ) 

Now they are:

global.i <- as.numeric(x%*%w%*%y) 
( global.p <- sum( abs(isim) > abs( global.i ) ) ) 

In one application I am doing, global.p is greater than 1.

Python version for spatialEco

Thank you for a great package! Are you planning on releasing a python version? (or did you?)

Alternatively, would you know of a python package computing spatial cross-correlations?

Thank you!

feature request? `only.holes`

Hi,

The remove.holes() functions is a handy utility function, but sometimes I need only the holes. A small change would all 'remove.holes()` to both remove or keep holes.

res <- lapply(1:length(xp), function(i) methods::slot(xp[[i]],  "Polygons")[!holes[[i]]])

to

res <- lapply(1:length(xp), function(i) methods::slot(xp[[i]],  "Polygons")[holes[[i]]])

Right now I just made my own only.holes function based on remove.holes.

Great package. Thanks!

Curvature Function

I was looking into the code for the curvature function, and I noticed a few issues:

1. For McNab and Bolstad, 36.2 m is hard-coded in as the distance to the outermost cell, but shouldn't this vary with cell size and focal window size?
2. For McNab, it seems that x[5] is hard-coded into the function. This would be the central cell in a 3x3 window, but not larger windows? Or am I missing something?

Thanks!

Bug in curvature function for mcnab type

Hi, I am helping a lab mate with some code and we're having issues with McNab's curvature metric. It seems to always return values equal to the input raster. Am I missing something?

library("terra")
library("spatialEco")

dem = terra::rast(system.file("extdata", "elev.tif", package = "spatialEco"))
mcnab = spatialEco::curvature(dem, type="mcnab")

terra::plot(dem-mcnab)

sessionInfo()

installation to get sample.distance

@jeffreyevans I have tried installing spatialEco (using remotes::install_github("jeffreyevans/spatialEco")) to get the version with sample. distance but I get the error sf is not found yet I have sf 1.0.6 version installed ( local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE) if(nchar(pkg)) library(pkg, character.only=TRUE)}) Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE Warning message: package ‘sf’ was built under R version 4.2.0 )

Kindly assist

`

remotes::install_github("jeffreyevans/spatialEco")
Downloading GitHub repo jeffreyevans/spatialEco@HEAD
Skipping 10 packages ahead of CRAN: numDeriv, colorspace, timeDate, zoo, sandwich, latticeExtra, abind, rgdal, multcomp, yaImpute

checking for file 'C:\Users\pc\AppData\Local\Temp\RtmpiKsqfA\remotes31d83a9d3d34\jeffreyevans-spatialEco-76d5739/DESCRIPTION' ...

checking for file 'C:\Users\pc\AppData\Local\Temp\RtmpiKsqfA\remotes31d83a9d3d34\jeffreyevans-spatialEco-76d5739/DESCRIPTION' ...

√ checking for file 'C:\Users\pc\AppData\Local\Temp\RtmpiKsqfA\remotes31d83a9d3d34\jeffreyevans-spatialEco-76d5739/DESCRIPTION'

• preparing 'spatialEco': (2.3s)
  checking DESCRIPTION meta-information ...

  checking DESCRIPTION meta-information ...

√ checking DESCRIPTION meta-information

• checking for LF line-endings in source and make files and shell scripts (377ms)

• checking for empty or unneeded directories

• building 'spatialEco_1.3-8.tar.gz'

Installing package into ‘C:/Users/pc/Documents/R/win-library/4.1’
(as ‘lib’ is unspecified)

• installing source package 'spatialEco' ...
  ** using staged installation
  ** R
  ** data
  *** moving datasets to lazyload DB
  ** inst
  ** byte-compile and prepare package for lazy loading
  ** help
  *** installing help indices
  converting help for package 'spatialEco'
  finding HTML links ... done
  DAYMET_tiles html
  annulus.matrix html
  ants html
  background html
  bearing.distance html
  breeding.density html
  built.index html
  cgls_urls html
  chae html
  chen html
  class.comparison html
  classBreaks html
  collinear html
  combine html
  concordance html
  conf.interval html
  convexHull html
  cor.data html
  correlogram html
  cross.tab html
  crossCorrelation html
  csi html
  curvature html
  dahi html
  date_seq html
  daymet.point html
  daymet.tiles html
  dispersion html
  dissection html
  divergence html
  download.daymet html
  download.hansen html
  download.prism html
  effect.size html
  elev html
  erase.point html
  explode html
  extract.vertices html
  focal.lmetrics html
  fuzzySum html
  gaussian.kernel html
  geo.buffer html
  group.pdf html
  hexagons html
  hli html
  hli.pt html
  hsp html
  hybrid.kmeans html
  idw.smoothing html
  impute.loess html
  insert html
  insert.values html
  is.empty html
  is.whole html
  kde.2D html
  kendall html
  kl.divergence html
  knn html
  land.metrics html
  local.min.max html
  loess.boot html
  loess.ci html
  logistic.regression html
  max_extent html
  moments html
  morans.plot html
  mwCorr html
  nni html
  nth.values html
  o.ring html
  oli.asw html
  optimal.k html
  optimized.sample.variance html
  outliers html
  overlap html
  parea.sample html
  parse.bits html
  partial.cor html
  plot.effect.size html
  plot.loess.boot html
  point.in.poly html
  poly.regression html
  polyPerimeter html
  polygon_extract html
  pp.subsample html
  print.cross.cor html
  print.effect.size html
  print.loess.boot html
  proximity.index html
  pseudo.absence html
  pu html
  quadrats html
  random.raster html
  raster.Zscore html
  raster.change html
  raster.deviation html
  raster.downscale html
  raster.entropy html
  raster.gaussian.smooth html
  raster.invert html
  raster.kendall html
  raster.mds html
  raster.modified.ttest html
  raster.moments html
  raster.transformation html
  raster.vol html
  rasterCorrelation html
  rasterDistance html
  remove.holes html
  rm.ext html
  rotate.polygon html
  sa.trans html
  sample.annulus html
  sample.line html
  sample.poly html
  sampleTransect html
  sar html
  se.news html
  separability html
  sg.smooth html
  shannons html
  shift html
  similarity html
  smooth.time.series html
  sobal html
  sp.kde html
  sp.na.omit html
  spatial.select html
  spectral.separability html
  spherical.sd html
  srr html
  stratified.random html
  subsample.distance html
  summary.cross.cor html
  summary.effect.size html
  summary.loess.boot html
  swvi html
  time_to_event html
  topo.distance html
  tpi html
  trasp html
  trend.line html
  tri html
  vrm html
  winsorize html
  wt.centroid html
  zonal.stats html
  *** copying figures
  ** building package indices
  ** testing if installed package can be loaded from temporary location
  *** arch - i386
  Error: package or namespace load failed for 'spatialEco' in library.dynam(lib, package, package.lib):
  DLL 'sf' not found: maybe not installed for this architecture?
  Error: loading failed
  Execution halted
  *** arch - x64
  ERROR: loading failed for 'i386'
• removing 'C:/Users/pc/Documents/R/win-library/4.1/spatialEco'
• restoring previous 'C:/Users/pc/Documents/R/win-library/4.1/spatialEco'
  Warning message:
  In i.p(...) :
  installation of package ‘C:/Users/pc/AppData/Local/Temp/RtmpiKsqfA/file31d867556dc4/spatialEco_1.3-8.tar.gz’ had non-zero exit status

`

proximity.index

Hi @jeffreyevans

I just started to use your package to calculate proximity.index but I only get zeros after calculations. Do I need to transform the CRS to UTM?

Any help would be greatly appreciated
Cheers,
Edicson

pol.cqf = st_read("poly_sepa/cqf1.shp")
mapview(pol.cqf, layer = "Admin. level 0", color = "#000000") #All good
px.cqf1 <-proximity.index(pol.cqf, min.dist = 3)
px.cqf1
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[47] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
:
:

My polygon's information

pol.cqf
class : SpatialPolygonsDataFrame
features : 463
extent : -74.36225, -74.264, 4.5605, 4.6605 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 2
names : DN, area
min values : 1, 766.881
max values : 1, 46445029.38

#Example polygon

meuse_poly
class : SpatialPolygonsDataFrame
features : 155
extent : 178562.3, 181429, 329677.5, 333650.5 (xmin, xmax, ymin, ymax)
crs : NA
variables : 13
names : cadmium, copper, lead, zinc, elev, dist, om, ffreq, soil, lime, landuse, dist.m, LU
min values : 0.2, 14, 37, 113, 5.18, 0, 1, 1, 1, 0, Aa, 10, forest
max values : 18.1, 128, 654, 1839, 10.52, 0.880389, 17, 3, 3, 1, W, 1000, nonforest

Incorporation of interaction effects in logistic.regression()

Hi, I am trying to fit a spatial binomial glm with interactions, but the arguments for thelogistic.regression() function only allow to specify the columns and not provide a formula. I assumed that I could center my covariates and manually create columns representing the interaction (e.g. taking the product of them) before adding them in as separate covariates.

My question is, is this the correct way to incorporate interactions, given we can't specify a formula?

The final code looks like this:

## dt is the table with the centered covariates. 
## `dist` is a binary response

## non-spatial glm
# mod <- glm(dist ~ cov1 + cov2 + cov3 + cov4 + 
#              cov5*cov6 + cov5*cov7 + cov5*cov8 +
#              cov5*cov9, family = binomial(link = "logit"), data=dt)

spo <- st_as_sf(dt, coords=c("x", "y"), crs=4326)
dt[, `:=` (cov56 = cov5 * cov6,
              cov57 = cov5 * cov7,
              cov58 = cov5 * cov8,
              cov59 = cov5 * cov9)]
out <- logistic.regression(dt, y="dist", x=c("cov1", "cov2", "cov3",  "cov4", "cov5", "cov6",
                                                "cov7", "cov8", "cov9", "cov56", "cov57", 
                                                "cov58", "cov59"),
                           autologistic=TRUE, coords=st_coordinates(spo), longlat=TRUE)

package ‘spatialEco’ is not available (for R version 3.6.3)

Hello everyone,
When I try to install the spatialEco package, I receive the below error in r:

package ‘spatialEco’ is not available (for R version 3.6.3)

Do you know how can I install it?
Thanks,

Error in breeding.density

When running breeding.density from the CRAN version I get the following error:

Error in Ops.units(di, bw) : 
  both operands of the expression should be "units" objects

I believe this is due to changes in the class of the object returned by sf::st_distance (sf version 1.0-9). I resolved it for now by changing the point.density function defined within breeding.density as follows:

point.density <- function(x, bw, self = TRUE) {
    d <- sf::st_distance(x, x)
    d <- matrix(as.numeric(d), nrow(d), ncol(d))

I can submit a pull request if you'd like, but I'm not sure this is actually the best solution. I installed from github and get the same error.

crossCorrelation - spatial weight matrix is not symmetric by default

The reference states the weight matrix is symmetric. (Chen pg.3) The function row-standardizes when scale.matrix is set to true, which violates the assumption. Possible fix is unity standardization as in w <- as.matrix(w/sum(w)).

st_geos_binop error in spatialEco::point.in.poly

Code :

track_data <- read_excel(paste(tracks_dir, tracks, sep = "/"))
check <- tracks
skipper <- file_path_sans_ext(tracks)

track_lats <- (track_data[1:nrow(track_data)- 1, "latitude"] + track_data[2:nrow(track_data), "latitude"]) / 2
track_lons <- (track_data[1:nrow(track_data)- 1, "longitude"] + track_data[2:nrow(track_data), "longitude"]) / 2

track_delta_ts <- track_data[1:nrow(track_data)- 1, "Time_s"]
track_speeds <- (track_data[1:nrow(track_data)- 1, "speed_kts"]) * 1.852 / 3.6

track_analysis_mk <- data.frame(track_lats, track_lons, track_delta_ts, track_speeds)
long <- track_analysis_mk$longitude
lat <- track_analysis_mk$latitude

sp::coordinates(track_analysis_mk) <- c("longitude", "latitude")
sp::proj4string(track_analysis_mk) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
track_analysis_mk <- sp::spTransform(track_analysis_mk, proj4string(SCANS_3))
track_analysis_mk_SCANS_3 <- spatialEco::point.in.poly(track_analysis_mk, SCANS_3, sp = TRUE, duplicate = TRUE)

Error :

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, (Simulation_VALS_2008-2021_for_PC.R#198): st_crs(x) == st_crs(y) is not TRUE
Show stack trace

PS:

track_analysis_mk is a SpatialPointsDataFrame
SCANS is a SpatialPolygonDataFrame

I run the code on a Window workstation and I have this error. If I run the same code on my macOS, the code runs well.

Function crossCorrelation - Change of code arrangement?

Hi!
First of all, thanks for the immediate change considering my first issue that I posted. Testing the change, i came across another issue that might be worthwhile solving.
Whilst testing the change, I provided the function with a prespecified weight matrix. I took me some time to realize that I had to set the dist.function argument to some arbitrary "none" to stop the function from altering the provided weight matrix.
As I was looking at the code, it seems as if this issue could be avoided since I think the dist.functions are only necessary if there is a new weight matrix to generate. So i changed the code arrangement to avoid an alteration of the prespecified weight matrix no matter the dist.function argument given. At least for the one example I tested it seemed to work fine.

This is the old arrangement (lines 109 - 136) of the line code

if( is.null(w) ) {
if( is.null(coords) ) stop("If no Wij matrix is provided, a coordinates matrix is required")
w <- sp::spDists( coords )
} else {
if(!class(w)[1] == "matrix") stop("Spatial weights must be in matrix form")
if(ncol(w) != length(x) | nrow(w) != length(x)) stop("Spatial weights matrix must be symmetrical and match x")
w[which(is.na(w))] <- 0
}
if( dist.function == "inv.power" ) {
message("Calculating spatial weights matrix using inverse power function")
w <- 1 / w
diag(w) <- 0
w <- w / sum(w)
} else if (dist.function == "neg.exponent") {
message("Calculating spatial weights matrix using negative exponent")
diag(w) <- NA
mu <- mean(w, na.rm=TRUE)
for(i in 1:nrow(w)) {
for(j in 1:nrow(w)) {
w[i,j] <- round(exp( (-2 * w[i,j]) / mu ),6)
}
}
diag(w) <- 0
} else if (dist.function == "none") {
message("Wij matrix is being left raw")
} else {
stop("Not a valid matrix option")
}

This would be the altered arrangement:

if( is.null(w) ) {
if( is.null(coords) ) stop("If no Wij matrix is provided, a coordinates matrix is required")
w <- sp::spDists( coords )
if( dist.function == "inv.power" ) {
message("Calculating spatial weights matrix using inverse power function")
w <- 1 / w
diag(w) <- 0
w <- w / sum(w)
} else if (dist.function == "neg.exponent") {
message("Calculating spatial weights matrix using negative exponent")
diag(w) <- NA
mu <- mean(w, na.rm=TRUE)
for(i in 1:nrow(w)) {
for(j in 1:nrow(w)) {
w[i,j] <- round(exp( (-2 * w[i,j]) / mu ),6)
}
}
diag(w) <- 0
} else if (dist.function == "none") {
message("Wij matrix is being left raw")
} else {
stop("Not a valid matrix option")
}
} else {
if(!class(w)[1] == "matrix") stop("Spatial weights must be in matrix form")
if(ncol(w) != length(x) | nrow(w) != length(x)) stop("Spatial weights matrix must be symmetrical and match x")
w[which(is.na(w))] <- 0
}

all raster commands are showing error in spatialEco package

I am not able to run data smoothing command : smooth.time.series(raster_stack), sg.smooth() and few others. It is showing error:
Error in .overlayList(x, fun = fun, filename = filename, forcefun = forcefun, :
cannot use this formula, probably because it is not vectorized

Please help. I tried to run the command with raster stack and brick both but non are working. And I tried package installation from gitHub and CRAN both. Still same error.

Add 2 independent variables

Hi,
Currently, I am working on spatial downscaling of LST data with the use of more than one explanatory variables. Could you tell me in what way I could implement more independent variables as 'x'?

Thank you,

Paulina

Please attend to retirement of rgdal, rgeos and maptools

See: https://r-spatial.org/r/2022/04/12/evolution.html. sp on CRAN now offers control of the environment variable, so the attached log is from _SP_EVOLUTION_STATUS_=2 _R_CHECK_FORCE_SUGGESTS_=FALSE R CMD check spatialEco_1.3-7.tar.gz:
00check.log

Unconditional use of package `rms`

spatialEco only Suggests the package rms, but R CMD check fails without this package installed since it is called unconditionally in the examples of logistic.regression. Maybe you would consider wrapping the example in a condition that ensures that it only runs if rms is present?

Bug when typing : library(spatialEco)

Dear Jeffrey,
I tried to install the library on my pc , I have the following problem ?
Any idea ?
Best, Mathieu Fallet

install.packages("spatialEco")
Installing package into ‘C:/Users/mathi/Documents/R/win-library/4.0’
(as ‘lib’ is unspecified)
essai de l'URL 'https://cran.rstudio.com/bin/windows/contrib/4.0/spatialEco_1.3-7.zip'
Content type 'application/zip' length 1193306 bytes (1.1 MB)
downloaded 1.1 MB

package ‘spatialEco’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\mathi\AppData\Local\Temp\RtmpGME2md\downloaded_packages

library(spatialEco)
Erreur : package or namespace load failed for ‘spatialEco’:
.onLoad a échoué dans loadNamespace() pour 'units', détails :
appel : udunits_init(path)
erreur : la fonction 'Rcpp_precious_remove' n'existe pas dans le package 'Rcpp'

polyPerimeter doesn't mention units

The function seems to work fine, I just have no idea what the units are that it is outputting to. My polygon output is in Km.sq.

An update on the help file would be perfect.

Best,

RG

Multiply a satellite image with gaussian kernel

I am trying to convolve a satellite image with gaussian kernel using the function spatialEco::gaussian.kernel but I am receiving this error: Error in gm %*% m : non-conformable arguments no matter the n value. Can you please help me?

Here is the code:

library(raster)
library(spatialEco)

pan = raster("C:/Users/Geography/Desktop/esda/mumbai/2014/nov/pan.tif")

m = as.matrix(pan)

gm = spatialEco::gaussian.kernel(sigma = 0.5, n = 5)

wm = gm %*% m

and here is the image:

new("RasterLayer", file = new(".RasterFile", name = "C:\\Users\\Geography\\Desktop\\esda\\mumbai\\2014\\nov\\pan.tif", 
    datanotation = "FLT4S", byteorder = "little", nodatavalue = -Inf, 
    NAchanged = FALSE, nbands = 1L, bandorder = "BIL", offset = 0L, 
    toptobottom = TRUE, blockrows = 256L, blockcols = 256L, driver = "gdal", 
    open = FALSE), data = new(".SingleLayerData", values = logical(0), 
    offset = 0, gain = 1, inmemory = FALSE, fromdisk = TRUE, 
    isfactor = FALSE, attributes = list(), haveminmax = TRUE, 
    min = 0.066125728189945, max = 0.30304148793221, band = 1L, 
    unit = "", names = "pan"), legend = new(".RasterLegend", 
    type = character(0), values = logical(0), color = logical(0), 
    names = logical(0), colortable = logical(0)), title = character(0), 
    extent = new("Extent", xmin = 582700, xmax = 604300, ymin = 1005700, 
        ymax = 1047400), rotated = FALSE, rotation = new(".Rotation", 
        geotrans = numeric(0), transfun = function () 
        NULL), ncols = 216L, nrows = 417L, crs = new("CRS", projargs = "+proj=lcc +lat_0=18.88015774 +lon_0=76.75 +lat_1=16.625 +lat_2=21.125 +x_0=1000000 +y_0=1000000 +datum=WGS84 +units=m +no_defs"), 
    history = list(), z = list())

I managed to solve the issue using Google Earth Engine but I would like to do the same in R.

sf support

Hello,

Can I use the raster.modifed.ttest() function for a regular gridded "sf" object?

bug in VRM?

I think the cos and sin in L33-34 are flipped from what they should be. These are the xy and z vectors.

currently, I get substantially different results from SAGA GIS but if I flip the definitions for the xy and z vectors then this function and SAGA implementaion are in agreement with each other.

spatialEco/R/vrm.R

Lines 33 to 34 in 8fc07e8

	sin.slp <- raster::calc(sa[["slope"]], fun=cos) # xyRaster
	cos.slp <- raster::calc(sa[["slope"]], fun=sin) # zRaster

Error between different packages with MannKendall test

Cheers!

I tested recently the MannKendall trend test. If on a limited dataset from one pixel (extracted on some point) I got the same results with 3 test from 3 different packages, when I try to compute the test on the whole rasters I got wrong values. Just the tau test value is the same with the other tests, and somehow the Pvalue is the one for the Z test value and the Intercept value the one for the Pvalue.
Anyone got any idea what can it be?

A transfer link for the raster test data: https://we.tl/t-FomZSMXStk

Test data for the for the different trend test applied in R.
test.csv

Screenshot from QGIS with all the values for the interest raster resulted after MannKendall test

Screenshot with the results from R and the example code with all the packages imported.

R snippet to apply MannKendall test from different packages

library(stringi)
library(Kendall)
library(dplyr)
library(trend)
library(EnvStats)



data22<- read.csv ("E:/stats_tile/test.csv", encoding = "UTF-8")##locale=locale(encoding = "UTF-8"))

test_1<-mk.test(data22$valori)
test_2<-MannKendall(data22$valori)
test_3<-spatialEco::kendall(data22$valori)
test_4<-EnvStats::kendallTrendTest(data22$valori)

test_1
test_2
test_3
test_4

tpi annulus option

Would it be possible to include an option for an annulus focal window in tpi.R?

Looks like Jan Holstein has attempted this in his crecs pacakge - however I'd prefer to use spatialEco since its on CRAN: https://rdrr.io/github/janhoo/crecs/src/R/terrain.R

Thank you for your attention to this.

class.comparison Raster error

Hello,

I am trying to use the class.comparison function to perform t.test on two landcover classification rasters in GeoTiff format. Now the ?class.comparison tells me that the function does support raster objects, but it's still giving me the error below. How can I fix this? Or do I have to convert the raster into a SpatialPixelsDataFrame?

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘gridded’ for signature ‘"RasterLayer"’

zonal.stats speed boost!

An issue and a suggestion all wrapped in one!

I've had the great fortune of finding spatialEco very useful for some zonal summaries I've been performing. However, it appeared that the rasterize() function was a bit of a processing speed bottleneck. Recently, a colleague of mine pointed me toward the fasterize library. It's designed to do basically one thing: make rasterizing vectors more efficient. I took the liberty of grabbing the zonal.stats() function code and came up with a slight edit, utilizing fasterize.

zonal.stats <- function(x, y, stat, trace = TRUE, plot = TRUE) {
  if (class(y) != "RasterLayer") 
    stop("y must be a raster object")
  if (class(x) != "SpatialPolygonsDataFrame") 
    stop("x must be a SpatialPolygonsDataFrame object")
  results <- vector()
  for (j in 1:nrow(x)) {
    if (trace == TRUE) {
      cat("processing", j, "of", nrow(x), "\n")
    }
    lsub <- sf::st_as_sf(x[j, ])
    cr <- raster::crop(y, raster::extent(lsub), snap = "out")
    crop.NA <- raster::setValues(cr, NA)
    fr <- fasterize::fasterize(lsub, cr)
    r <- raster::mask(x = cr, mask = fr)
    if (plot == TRUE) {
      plot(r, main = paste("Polygon: ", j, sep = " "))
    }
    r <- raster::values(r)
    r <- stats::na.omit(r)
    if (length(r) < 1) {
      results <- append(results, NA)
    } else {
      results <- append(results, stat(r))
    }
  }
  results
} 

In my testing (400 NPS polygons, 90 m raster) this change dropped the processing time from ~ 60 min to ~ 3 min. That's the only dataset I've tested it on but I'd assume similar results would be evident in other examples.

In any case, I figured I would drop this here to see if you felt it was appropriate to incorporate.

Cheers!

Peder

Error: Loading required package: spatialEco

Hello

Installation of spatialEco in a docker container running ubuntu succeeds. I also install sf, regeos, spdep and terra just in case.

RUN R -e "install.packages('sf')"
RUN R -e "install.packages('regeos')"
RUN R -e "install.packages('spdep')"
RUN R -e "install.packages('terra')"
RUN R -e "install.packages('spatialEco')"

But when it comes to loading and using the package, it fails. Here is the exact error

Loading required package: spatialEco
Error: package or namespace load failed for ‘spatialEco’ in dyn.load(file, DLLpath = DLLpath, ...):
	unable to load shared object '/usr/local/lib/R/site-library/terra/libs/terra.so':
libproj.so.15: cannot open shared object file: No such file or directory
Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘listenv’, ‘parallelly’, ‘future’, ‘globals’, ‘future.apply’, ‘progressr’, ‘numDeriv’, ‘SQUAREM’, ‘lava’, ‘colorspace’, ‘prodlim’, ‘utf8’, ‘bit’, ‘prettyunits’, ‘png’, ‘jpeg’, ‘RColorBrewer’, ‘viridisLite’, ‘farver’, ‘labeling’, ‘munsell’, ‘matrixStats’, ‘evaluate’, ‘highr’, ‘yaml’, ‘xfun’, ‘backports’, ‘gower’, ‘ipred’, ‘lubridate’, ‘tidyr’, ‘timeDate’, ‘ellipsis’, ‘pkgconfig’, ‘vctrs’, ‘fansi’, ‘pillar’, ‘bit64’, ‘tidyselect’, ‘withr’, ‘progress’, ‘Formula’, ‘latticeExtra’, ‘gtable’, ‘gridExtra’, ‘data.table’, ‘viridis’, ‘base64enc’, ‘isoband’, ‘scales’, ‘MatrixModels’, ‘conquer’, ‘mvtnorm’, ‘TH.data’, ‘sandwich’, ‘knitr’, ‘checkmate’, ‘htmlwidgets’, ‘rstudioapi’, ‘fastmap’, ‘xts’, ‘intervals’, ‘ModelMetrics’, ‘pROC’, ‘recipes’, ‘iterators’, ‘generics’, ‘gibble’, ‘purrr’, ‘decido’, ‘unjoin’, ‘gridBase’, ‘crsmeta’, ‘cli’, ‘clipr’, ‘hms’, ‘tibble’, ‘vroom’, ‘tzdb’, ‘bitops’, ‘Hmisc’, ‘ggplot2’, ‘SparseM’, ‘quantreg’, ‘polspline’, ‘multcomp’, ‘htmlTable’, ‘htmltools’, ‘digest’, ‘fastmatrix’, ‘nortest’, ‘zoo’, ‘spacetime’, ‘FNN’, ‘caret’, ‘XML’, ‘geosphere’, ‘rgdal’, ‘doParallel’, ‘foreach’, ‘dplyr’, ‘fasterize’, ‘silicate’, ‘exactextractr’, ‘readr’, ‘RCurl’, ‘RANN’, ‘rms’, ‘yaImpute’, ‘SpatialPack’, ‘EnvStats’, ‘maptools’, ‘GeNetIt’, ‘gstat’, ‘RStoolbox’, ‘tabularaster’

What do I do to get it to work? R version is 4.0.5.
Any ideas. Help please.

x and y comment flipped in VRM

This does not affect results, but in the VRM function, the lines commented as xraster and yraster are switched. sin.asp should be "#xraster" and cos.asp should be "#yraster". They are used correctly in subsequent lines that use x.sum and y.sum (i.e. x.sum is calculated using sin.asp which is xrast even though it is commented as yrast).

Comparing (gaussian) smoothing algorithms - spatialEco vs terra

I am comparing smoothing algorithms in R. More specifically, the raster.gaussian.smooth from the spatialEco package and the focal from the terra package. Visually, the results look similar but the pixel values are way different. It seems that raster.gaussian.smooth rescales the data or something. In the picture below, the first layer is my original raster, the second is the result after using raster.gaussian.smooth and the third is the result from terra's focal:

How can I stop raster.gaussian.smooth from rescaling (or whatever it does) the data? is this possible?

Here is the code:

library(terra)
library(spatialEco)

r = rast("path/tirs.tif")

# spatialEco
smoothed = spatialEco::raster.gaussian.smooth(r, sigma = 0.5 * 350)

# terra focal
gf <- focalMat(r, 0.5*350, "Gauss")
r_gf <- focal(r, w = gf)

Here is the raster I am using:

tirs = rast(ncols=578, nrows=449, nlyrs=1, xmin=504400, xmax=562200, ymin=155700, ymax=200600, names=c('B10_median'), crs='EPSG:27700')

asc.from.raster gone since v1.3-0

Hi folks,

I'm looking to run some old code from a friend, which uses asc.from.raster. This function appeared in spatialEco v1.3-0, and looks to still be in SDMTools, however their latest version isn't available for the latest R.

Do you remember the rationale for removing asc.from.raster?

I don't suppose you know where I could find that code? I was digging through the commits but feels like I could be here a while!

Cheers!

Simon

focal.lmetrics output doesn't include 0 or 1

I'm trying to output a grid that captures that proportion of forest (value 1) in a landscape. The code below works to accomplish that; however, where the % forest is 0% or 100%, no data is written to the output grid. Is there anyway to output these values to the output grid like the moving window feature of FRAGSTATS does? Or do I have to use FRAGSTATS for this? Thanks.

land_forest <- raster("gis/moving_window/land_use_grid_clips/14.img")
land_forest_moving <- focal.lmetrics(land_forest, w = 11, land.value = 1,
                                       metric = "prop.landscape")
writeRaster(land_forest_moving, "gis/moving_window/outputs/14_forest.img")

rasterCorrelation producing rings in output raster

I am reading in two .tif files as r1 and r2 using raster::raster("*.tif") and trying to do spatial correlation between the two by using the following command.
r.cor <- rasterCorrelation(r1, r2, s = 5, type = "pearson")
First, the program is extremely slow and takes about 20 minutes for a relatively small raster. Second, I am getting circular rings in my output raster. The example provided with the function do not produce any of such affect.

wt.centroid and tibbles

I have just answered this question on stackoverflow https://stackoverflow.com/questions/76334169/how-to-compute-the-weighted-centroid-of-multiple-st-points-with-the-spatialeco-p .

There seems to be a problem in the code for wt.centroid as the method it uses to extract the column of weights does not work if the input sf object is also a tibble. It seems that sf::st_drop_geometry behaves differently for dataframes and tibbles (so perhaps I ought to be raising the issue there).

Regarding negative p-value for raster.kendall

Thank you for this great package. I was trying to run Mann Kendall test on the raster stack. While I was looking at the p-value of the output. I am getting negative values and many values greater than 1. I believe there is some issue in it. It would be great if you look into that.

Thanks

raster.kendall function in spatialEco package

I have 6 climate layer and stack them as a file. when I run "raster.kendall" function, I get an Error:
"Error in raster.kendall(x, p.value = TRUE, z.value = TRUE, intercept = TRUE, :
please install EnvStats package before running this function"

When I install "EnvStat" package, and run again "raster. kendall" function, I get another error:
"Error in (function (y, x = seq(along = y), alternative = "two.sided", : When ci.slope=TRUE, there must be at least 3 non-missing, finite observations"

crossCorrelation - sum(lsci.xy) != sum(lsci.yx)

The sums of lsci vectors are not equal to each other.
Gsci value is divided by n-1.

# data
data(meuse)
x = meuse$zinc
y = meuse$copper
x = (x - mean(x))/(stats::sd(x) * sqrt((length(x) - 1)/length(x)))
y = (y - mean(y))/(stats::sd(y) * sqrt((length(y) - 1)/length(y)))

# LSCI calculation
w = as.matrix(dist(cbind(meuse$x, meuse$y)))
w = 1/(1+w)
diag(w) = 0
w = w/sum(w)
l.xy = x * t(y) %*% w
l.yx = y * t(x) %*% w

#print GSCI
sum(l.xy) # gsci
sum(l.yx) # same as above

# LSCI calcalation by spatialEco
cross.cor.spatialEco = crossCorrelation(x = x, y = y, w = w, scale.matrix = FALSE)
cross.cor.spatialEco$I # equal to sum(l.xy)/(n-1)
sum(cross.cor.spatialEco$SCI$lsci.xy) 
sum(cross.cor.spatialEco$SCI$lsci.yx) # different from above

Trasp function

Hi there,
I think that there is missing parenthesis in the equation.
ORG:
"trasp <- function(x, ...) {
if (!inherits(x, "RasterLayer")) stop("MUST BE RasterLayer OBJECT")
asp <- raster::terrain(x, out='aspect', unit='degrees')
return( raster::calc(asp, fun=function(x) { 1 - cos( (3.142/180) *(x - 30) ) / 2 }, ... ) )"

Should not be?:

"trasp <- function(x, ...) {
if (!inherits(x, "RasterLayer")) stop("MUST BE RasterLayer OBJECT")
asp <- raster::terrain(x, out='aspect', unit='degrees')
return( raster::calc(asp, fun=function(x) { (1 - cos( (3.142/180) *(x - 30) )) / 2 }, ... ) )"

Best,
Monika

Include sf and rgdal as dependency

Hi @jeffreyevans

Thank you for making this package, really helpful. I'm currently using it for calculating VRM. But there's a minor issue that you may want to address.

1. Include sf when installing this package so that other users don't have to install it manually if that's possible

> library(raster)
Loading required package: sp
> library(spatialEco)
Error: package or namespace load failed for ‘spatialEco’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
 there is no package called ‘sf’

2. Include rgdal when installing this package so that other users don't have to install it manually if that's possible

> data(elev)
> elev <- projectRaster(elev, crs="+proj=robin +datum=WGS84",res=1000, method='bilinear')
Error in .requireRgdal() : package 'rgdal' is not available

Because in this case I need to calculate vrm which needs elev as an example. This works out after all, but automatically installing the dependency may help other users perhaps.

Thanks once again :)

sf.kde flips raster output

Dear Jeff,

The sf.kde() functions gives me a raster output with a different orientation. The output does not match the locations of the input vector features. Here is my code.

library(spatstats)
library(terra)

data("nbfires")

# convert ppp to sf
nbfires <- rescale(nbfires)
nbfires <- subset(nbfires, fnl.size > 0)
boundary <- st_as_sf(nbfires$window)
df.fires <- data.frame(x=nbfires$x, y=nbfires$y, marks(nbfires))
fires <- st_as_sf(df.fires, coords=c("x", "y"))

fire_count <- terra::rasterize(fires, ras, fun = "length", background = 0)
plot(fire_count)

fire_kde <- spatialEco::sf.kde(fires, bw = c(50, 50), ref = fire_count, standardize = T)

# wrong orientation
plot(fire_kde)
plot(boundary, add=T)

# right orientation - but does not match
plot(flip(flip(t(fire_kde), 'horizontal'), 'vertical'))
plot(boundary, add=T)

Thank you!
Dirk