R: Вычислить большие расстояния вдоль меридианов в метрическом формате. - PullRequest
Новая
       37

R: Вычислить большие расстояния вдоль меридианов в метрическом формате.

0 голосов
/ 27 марта 2019

Я хотел бы вычислить расстояния от всех пикселей в глобальном растре с 5-минутным разрешением до экватора.Полученные расстояния должны быть в метрах с учетом кривизны Земли и не должны искажаться проекциями карты.

Приложения на небольших территориях (например, в отдельных странах) обычно преодолевают угрозу искажений, вызванных проекцией карты, путемвыбрать код EPSG, который хорошо подходит для данного региона.Поскольку мое приложение охватывает весь земной шар, этот трюк не работает.Однако мне не нужна карта, которая адекватно отображает планету во всех аспектах.Мне нужна только проекция, которая не искажает расстояние между каждым пикселем и ближайшей точкой на экваторе, т.е. проекция, которая не искажает расстояния вдоль меридианов.Насколько мне известно, проекция Плато Карре ("+proj=longlat +datum=WGS84") удовлетворяет этому критерию.Учитывая эту информацию, я попробовал следующий код: 

r_res <- 1/12
r <- raster(resolution = c(r_res, r_res), crs = "+proj=longlat +datum=WGS84")
p <- as(r, "SpatialPoints")
equator <- st_sfc(st_point(c(-180,0)), st_point(c(0,0)), st_point(c(180,0))) %>% st_combine() %>% st_cast(., "LINESTRING") %>% st_sf(., crs = 4326) %>% st_transform(crs = "+proj=longlat +datum=WGS84") %>% as(., "Spatial")
d <- gDistance(p, equator, byid = T)
dmin <- apply(d, 2, min)
r[] <- dmin

Расстояния, вычисленные в этом примере, к сожалению, выражены в градусах, а не в метрах, что вызвано форматом проекции longlat.Кроме того, я не смог найти никакой информации о том, учитывает ли rgeos::gDistance() кривизну планеты.И согласно соответствующей публикации, gDistance() не следует даже применять к longlat данным.

В качестве альтернативы, я спроецировал сетку в другую проекцию (Mollweide), которая производит вывод в метрах: 

r_res <- 1/12
r <- raster(resolution = c(r_res, r_res), crs = "+proj=longlat +datum=WGS84")
r <- projectRaster(r, crs="+proj=moll +ellps=WGS84")
p <- as(r, "SpatialPoints")
equator <- st_sfc(st_point(c(-180,0)), st_point(c(0,0)), st_point(c(180,0))) %>% st_combine() %>% st_cast(., "LINESTRING") %>% st_sf(., crs = 4326) %>% st_transform(crs="+proj=moll +ellps=WGS84") %>% as(., "Spatial")
d <- gDistance(p, equator, byid = T)
dmin <- apply(d, 2, min)
r[] <- dmin

В этом случае, однако, я не уверен, исправляет ли gDistance() искажения расстояния, подразумеваемые проекцией Молвейде.

Мне известно, что расстояния от экватора могут быть рассчитаны с использованием правилабольшой палец: latitude * 111 km.Это приближение, однако, слишком неточно для моих приложений, которые требуют функции более высокой точности.

Было бы здорово, если бы кто-нибудь мог дать какой-то совет.Не стесняйтесь полагаться на функции расстояния, отличные от gDistance().Пока расстояния не искажены, измеряются в метрах и учитывают кривизну земли, я в порядке с любой функцией, будь то от sf, gdistance, raster, geosphere, rgeosили любой другой пакет.

1 Ответ

2 голосов
/ 27 марта 2019

Вот подход (который не безопасен для памяти для очень больших растров) 

library(raster)
r_res <- 1/12
r <- raster(resolution = c(r_res, r_res), crs = "+proj=longlat +datum=WGS84")
# latitudes for one column
lats <- cbind(0, yFromRow(r, 1:nrow(r)))
#distances (in km) for one column
dist <- pointDistance(lats, cbind(0,0), lonlat=TRUE) / 1000
# assign to all cells
d <- setValues(r, rep(dist, each=ncol(r)))

Это установило бы "правило большого пальца" на mean(dist/abs(lats[,2])) = 110.8032 км на градус широты

Вот безопасный для памяти (но неэффективный) метод: 

x <- init(r, "y")
f <- function(i) { pointDistance(cbind(0, i), cbind(0,0), lonlat=TRUE) / 1000} 
z <- calc(x, fun=f)
...