Search code examples
rdistancerasterspatial

raster::distance() not recognizing NA cells

I am looking to calculate distances from NA cells to a spatial object. However, the distance() function is not recognizing NA cells that I have created. Here is reproducible code:

To make spatial object:

library(raster)
CAN <- getData("GADM", country = "Canada", level = 1)#level 1 just gives borders for provinces

library(sf)#to make sf object
sf_CAN <- as(CAN,"sf")

st_crop_CAN <- st_crop(sf_CAN, xmin = -60.17, xmax = -59.65, ymin = 43.92, ymax = 44.05) #bounding box gives sable island

Create raster file and use mask to then use distance():

library(marmap)
shelf <- getNOAA.bathy(lon1 = -75, lon2 = -45,
                       lat1 = 35, lat2 = 60, resolution = 4)
library(raster)
CAN <- getData("GADM", country = "Canada", level = 1)#level 1 just gives borders for provinces
library(sf)#to make sf object
sf_CAN <- as(CAN,"sf")

st_crop_CAN <- st_crop(sf_CAN, xmin = -60.17, xmax = -59.65, ymin = 43.92, ymax = 44.05) #bounding box gives sable island

shelf_dis <- marmap::as.raster(shelf)
values(shelf_dis) <- NA
newR <- mask(shelf_dis,as(st_crop_CAN$geometry, "Spatial"))
dist_sable <- distance(newR)

#error code after distance()

Error in .local(x, y, ...) : 
  RasterLayer has no NA cells (for which to compute a distance)

One issue that I see is that in the min and max values of newR object are NA as shown below:

> newR
class      : RasterLayer 
dimensions : 375, 450, 168750  (nrow, ncol, ncell)
resolution : 0.06666667, 0.06666667  (x, y)
extent     : -75, -45, 35, 60  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : NA, NA  (min, max)

So it is unable to distinguish between NA cells and the spatial object because there are only NA cells to begin with?

I am not sure if my approach to this is the best way. Any advice or help appreciated. At the end of the day I am just trying to get distances from points (not provided) to the spatial object. That part seems easy enough with extract(). I just need to get this new raster file! This is my first go at spatial work in R.

Solution

  • Yes, there must be cells that are not NA to compute distance to. The general approach is usually like this

    library(raster)
# example polygons
filename <- system.file("external/lux.shp", package="raster")
p <- shapefile(filename)   
pp <- crop(p, extent(5.698, 6.231, 49.442, 49.781))
# create or use existing raster
r <- raster(p, res=0.01)

r <- rasterize(pp, r)
d <- distance(r)

d
#class      : RasterLayer 
#dimensions : 73, 78, 5694  (nrow, ncol, ncell)
#resolution : 0.01, 0.01  (x, y)
#extent     : 5.74414, 6.52414, 49.45162, 50.18162  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : 0, 49114.11  (min, max)