I'm trying to do some raster algebra and one of my layers is a distance from a coastline layer. I'm using rasters with 12.5m resolution and when I use the distance function, I seem to run into memory issues. Is there a way to set a maximum distance like dist in gDistance? Are there other ways to speed this calculation up? Here's a mock raster with similar specs:

#build blank raster modeled after large raster
r <- raster(xmn=1035792, xmx= 1116792, ymn=825303.6, ymx=937803.6, resolution = 12.5,crs = "+init=epsg:3174")
r <- setValues(r, 0)

#create lake polygon
x <- c(1199999, 1080000, 1093067, 1090190, 1087977, 1070419, 1180419)
y <- c(957803.6,937803.6, 894366.9, 872153.9, 853703.0, 825353.6, 805353.6)
poly.lake <- SpatialPolygons(list(Polygons(list(Polygon(data.frame(x,y))), ID = 1)))

# make values NA where lake polygon does not intersect raster
r <- mask(r, poly.lake)

#run distance function
r.dist <-distance(r)
  • 1
    I wouldn't call 12.5 m "high resolution".
    – Iris
    Sep 14 '16 at 14:22
  • Is the point you're measuring to always a line? There may be a way to do that faster. Moreover, what kind of machine are you running? I just ran your code, took 30 minutes to complete, and only used 2Gb of RAM, if your system has less than 4 Gb of RAM, do you have access to a larger system? If your grids double, or triple in size you may need to consider chunking the work, which can get messy. Let me know what your specs are and how your data is typically arranged and I'll keep thinking. 12.5 m can be high resolution, if you're looking from the perspective of Mars. All about perspective.
    – Badger
    Sep 15 '16 at 20:54
  • @Badger, yes the dummy raster represents a lake polygon, which could be a line. I've tried the doEdge=TRUE argument, but it didn't help. My machine has 16 GB of RAM.
    – alaybourn
    Sep 16 '16 at 21:17
  • I tried running this again and it took ~ 70 minutes, but at least it ran. I'd still be interested in methods for quickening this calculation.
    – alaybourn
    Sep 19 '16 at 13:03

If you a using a lake boundary polygon in your analysis, your example should have something like that too.

Here is how I am interpreting your question:

I am creating a dummy lake polygon then masking it from a raster of zeros of the same specs you are using. Running distance on this set only took my machine 9 minutes.

You can reduce your search area by buffering the lake polygon/polygons, then use the buffer extent to clip your raster.

If you are interested in the interior distance of the lake to shore, use the inverse argument for the mask function. A code example is below.


# build blank raster - are you familar with the "+init=epsg:xxxx" string for projections?
r <- raster(xmn=1035792, xmx= 1116792, ymn=825303.6, ymx=937803.6, resolution = 12.5,crs = "+proj=aea +lat_1=42.122774 +lat_2=49.01518 +lat_0=45.568977 +lon_0=-84.455955 +x_0=1000000 +y_0=1000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
r <- setValues(r, 0)

# build lake polygon
x <- c(1088419, 1093067, 1090190, 1087977, 1088419)
y <- c(895030.8, 894366.9, 892153.9, 893703.0, 895030.8)
poly <- SpatialPolygons(list(Polygons(list(Polygon(data.frame(x,y))), ID = 1)))

# make values NA where polygon intesects raster
r <- mask(r, poly)

# run distance check
rD <- distance(r)

# to restrict your search
searchLimit <- 2500 # the maximum distance in raster units from lakes/boundaries
polyB <- gBuffer(poly, width = searchLimit, byid = T)
rB <- crop(r, extent(polyB))

# much faster...
rDB <- distance(rB)

Coming in late in the game, but try out the wbt_euclidian_distance function from the whitebox tools package in r. I had a large raster (12k x 8k resolution) that would choke on the distance function from the raster or terra packages. On the whitebox tools it took less than a second. Just remember that the wbt functions don't work on raster objects, but take input strings as paths to your raster object to be processed. Also, instead of coding your raster as 1s and NAs (where the raster::distance calculates the distance of NA cells to nonNA cells), you must code your raster values 1s and 0s.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.