R - найти самую дальнюю точку от множества точек на растеризованной карте США

Новое в пространственном анализе на R здесь. У меня есть шейп-файл для США, который я скачал с ЗДЕСЬ. У меня также есть набор точек широты и долготы (полмиллиона), которые лежат на территории США.

Я хотел бы найти «самую удаленную точку» — точку в пределах прилегающих США, которая дальше всего от набора точек.

Я использую пакеты rgdal, raster и sp. Вот воспроизводимый пример со случайной выборкой из 10 точек:

# 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)

После последних двух строк кода я получаю этот график.

построение расстояний от точек

Однако я бы хотел, чтобы это было поверх растровой карты США. И я бы хотел, чтобы он учитывал расстояния только от ячеек, которые находятся в смежных США, а не от всех ячеек в ограничивающей рамке. Как мне это сделать?

Я также был бы признателен за любые другие советы относительно файла формы, который я использую - лучший ли он? Должен ли я беспокоиться об использовании правильной проекции, поскольку мой фактический набор данных широта/долгота? Сможет ли DistanceFromPoints эффективно обрабатывать такой большой набор данных или есть функция получше?


person hoggue    schedule 07.07.2016    source источник


Ответы (2)


Чтобы ограничить растр D смежной территорией США, вы можете найти элементам rr присвоенные значения NA (т. е. растровые ячейки внутри ограничивающей рамки, но за пределами полигонов США) и присвоить этим же элементам D значение NA.

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

plot(D)
lines(usa)

Вы можете использовать «proj4string (usa)», чтобы найти информацию о проекции для шейп-файла usa. Если интересующие вас координаты основаны на другой проекции, вы можете преобразовать их, чтобы они соответствовали проекции шейп-файла США, следующим образом:

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

Не уверен в относительной эффективности DistanceFromPoints, но на моем компьютере для запуска на вашем примере с 10 точками потребовалось всего ~ 1 секунда.

person pbee    schedule 07.07.2016
comment
я не могу найти в документации {raster}, являются ли расстояния от DistanceFromPoints евклидовыми расстояниями или расстояниями по большому кругу. я предполагаю/надеюсь, что это последнее. он говорит, что для непроецированных координат долготы/широты это будет в метрах, но это все. ты знаешь? - person hoggue; 08.07.2016
comment
Пришлось заглянуть в источник {растра}... distFromPoints вызывает функцию distanceToNearestPoint, которую вы можете найти в /raster/src/distance.c. Если растр/точки не проецируются (широта), distanceToNearestPoint вызывает функцию distCos, которая возвращает расстояние по большому кругу на основе радиуса земли в метрах (источник в /raster/src/dist_util.c). - person pbee; 08.07.2016

Я думаю, вы искали функцию mask.

library(raster)
usa <- getData('GADM', country='USA', level=1)

# exclude Alaska and Hawaii
usa <- usa[!usa$NAME_1 %in% c( "Alaska" , "Hawaii"), ]

# get the extent and create raster with preferred resolution
r <- raster(floor(extent(usa)), res=1)
# rasterize polygons
rr <- rasterize(usa, r)

set.seed(89)
sample <- spsample(usa, 10, type = "random")

# Find distance from sample points to cells of USA raster
D <- distanceFromPoints(object = rr, xy = sample) 
# remove areas outside of polygons
Dm <- mask(D, rr)
# an alternative would be mask(D, usa)

# cell with highest value
mxd <- which.max(Dm)

# coordinates of that cell
pt <- xyFromCell(r, mxd)

plot(Dm)
points(pt)

Расстояния должны быть в порядке, также при использовании данных долготы/широты. Но rasterFromPoints действительно может быть немного медленным с большим набором данных, поскольку он использует алгоритм грубой силы.

person Robert Hijmans    schedule 08.07.2016