Я хотел бы вычислить расстояния от всех пикселей в глобальном растре с 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
или любой другой пакет.