EWS on subset data about spatialwarnings HOT 3 CLOSED

Hello and thanks for your interest in our package! I've given it some thoughts over the past few days.

Subsetting the raster data by itself is pretty standard (e.g. see ?raster::crop and ?raster::mask for examples), but computing EWS on nonrectangular data is a bit tricky right now with spatialwarnings. It is definitely something we should look into in the future though. We somehow need to find a way to tell spatialwarnings that missing values should be ignored for some computations, while still providing a sensible result.

At the moment many of the low-level functions that compute indicators in spatialwarningscannot handle NAs. For some of them as far as I know it is an intrinsic limitation (e.g. spectrum-based indicators), but for others it is just a technical limitation. For example there is no reason why we couldn't compute spatial variance while just ignoring NAs in the input matrix.

This can be done by telling spatialwarnings to use custom functions to compute indicators and test significance. I am giving you below bits of code you can adapt to your particular case. Please note that I pushed some changes in the package for this code to work (spatialwarnings in the current version will throw an error if it finds NAs in the original matrix), so you will need to install the development version.

I am assuming in what follows that your input data is a matrix with continuous values, where NAs represent values outside of the areas of interest (e.g. land). Here I am building such matrix for the example:

# You need the development version of spatialwarnings for this to work
library(devtools)
install_github("spatial-ews/spatialwarnings")
library(spatialwarnings)
#> Loading required package: future
#> This is spatialwarnings 2.99.99 (development version, use at your own risk!)
#> Use plan(multiprocess) or plan(multicore) to set up parallel processing

# Build and display an example matrix
example_mat <- outer(seq.int(100), seq.int(100), function(x, y) { 
  ifelse((x-40)^2 + (y-40)^2 < 10^2, TRUE, NA)
}) 
example_mat <- example_mat * runif(100*100)
display_matrix(example_mat)

We then define a custom 'indicator function' to compute the generic EWS on a matrix with missing values. We can call it on the original data using compute_indicator() (see also https://alex.lecairn.org/spatialwarnings-faq.html#how_can_i_add_another_indicator and ?compute_indicator for more information on custom indicator functions).

# This is our example function. cg_subsize is the coarse_grain length, see 
# ?coarse_grain for details
na_aware_ews <- function(mat, cg_subsize) { 
  require(moments) # required packages
  require(raster)
  mat_coarse <- coarse_grain(mat, cg_subsize)
  c(skewness = skewness(as.vector(mat_coarse), na.rm = TRUE), 
    variance = var(as.vector(mat_coarse), na.rm = TRUE), 
    moran    = Moran(raster(mat_coarse)))
}

indics <- compute_indicator(example_mat, fun = na_aware_ews, 
                            cg_subsize = 2)
#> Warning in check_mat(mat): NAs in provided matrix.
#> Loading required package: moments
#> Loading required package: raster
#> Loading required package: sp
#> 
#> Attaching package: 'raster'
#> The following object is masked from 'package:future':
#> 
#>     values
summary(indics)
#> Spatial Early-Warning: na_aware_ews 
#> 
#>  1 matrix (size: 100x100)
#> 
#>  Mat. # skewness variance moran
#>       1     0.18    0.021  0.03
#> 
#> The following methods are available:
#> as.data.frame indictest plot summary 

A tricky point arises if you want to test the significance of indicator values against randomized matrices. We need to make sure to produce null matrices that take into account the empty areas, i.e. keep NAs where NAs occur in the original matrix, and randomize the rest. The following function does just that. We pass it to indictest() to indicate we want to use it instead of the default canned methods.

# This is the function that given a matrix, returns a randomized matrix
randomize_matrix_no_na <- function(mat) { 
  mat[!is.na(mat)] <- sample(mat[!is.na(mat)])
  return(mat)
}

# Here our randomizing function keeps the NAs in place, instead of spreading 
# them all over the matrix. 
display_matrix(randomize_matrix_no_na(example_mat))

# Test significance
indics_test <- indictest(indics, nulln = 49, 
                         null_method = randomize_matrix_no_na)
summary(indics_test)
#> Spatial Early-Warning: na_aware_ews 
#> 
#>  1 matrix (size: 100x100)
#> 
#>  Mat # skewness P>null     variance P>null     moran P>null    
#>      1     0.16   0.28        0.018    0.6      0.12   0.02 *  
#> 
#>  Null model: custom function (49 null matrices used)
#>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> The following methods are available:
#> 

Note that in the above code you can substitute example_mat by a list of matrices (or raster objects). This will allow you to plot trends using plot().

I hope that helps with your problem, with the code above you should be able to compute generic indicators for each of the four areas, but of course it would have to be extended if you want to work with other indicators.

from spatialwarnings.

Hi there, When we were analysing data for Eby et al, 2017, GEB (https://doi.org/10.1111/geb.12570) - we faced similar problems with NA's. What we found that NA's were a tiny percentage of all data points (<5% or much less I think). Our data was binary -- 0 or 1 at each pixel/spatial location matrix. So we just replaced all NA's randomly with 1 (or 0) a probability equal to the proportion of 1's (or proportions of 0's) in the landscape. That kept the proportion of 1 (or 0) on landscape same as before.

One problem that may happen is if there are spatial correlations in the distribution of NAs. otherwise, the above simple method can be used as a pre-processing step before

from spatialwarnings.

Closing the discussion now.

from spatialwarnings.