R - find point farthest from set of points on rasterized USA map

Question

New to spatial analysis on R here. I have a shapefile for the USA that I downloaded from HERE. I also have a set of lat/long points (half a million) that lie within the contiguous USA.

I'd like to find the "most remote spot" -- the spot within the contiguous USA that's farthest from the set of points.

I'm using the rgdal, raster and sp packages. Here's a reproducible example with a random sample of 10 points:

# Set wd to the folder tl_2010_us_state_10
usa <- readOGR(dsn = ".", layer = "tl_2010_us_state10")

# Sample 10 points in USA
sample <- spsample(usa, 10, type = "random")

# Set extent for contiguous united states
ext <- extent(-124.848974, -66.885444, 24.396308, 49.384358)

# Rasterize USA
r <- raster(ext, nrow = 500, ncol = 500)
rr <- rasterize(usa, r)

# Find distance from sample points to cells of USA raster
D <- distanceFromPoints(object = rr, xy = sample) 

# Plot distances and points
plot(D)
points(sample)

After the last two lines of code, I get this plot.

However, I'd like it to be over the rasterized map of the USA. And, I'd like it to only consider distances from cells that are in the contiguous USA, not all cells in the bounding box. How do I go about doing this?

I'd also appreciate any other tips regarding the shape file I'm using -- is it the best one? Should I be worried about using the right projection, since my actual dataset is lat/long? Will distanceFromPoints be able to efficiently process such a large dataset, or is there a better function?

Answer 1

To limit raster D to the contiguous USA you could find the elements of rr assigned values of NA (i.e. raster cells within the bounding box but outside of the usa polygons), and assign these same elements of D a value of NA.

D[which(is.na(rr[]))] <- NA

plot(D)
lines(usa)

You can use 'proj4string(usa)' to find the projection info for the usa shapefile. If your coordinates of interest are based on a different projection, you can transform them to match the usa shapefile projection as follows:

my_coords_xform <- spTransform(my_coords, CRS(proj4string(usa)))

Not sure about the relative efficiency of distanceFromPoints, but it only took ~ 1 sec to run on my computer using your example with 10 points.