Широта и долгота, x и y, а также различия при их построении: геосфера и ggmap.

Я использую данные о широте и долготе и пытаюсь графически отобразить метрики в каждой точке карты Стокгольма (в зависимости от близости к этой точке). Меня больше интересует, чтобы точки были равномерно разнесены на изображении, а не одинаково разнесены по фактическому расстоянию: в этом смысле я понимаю, что расстояния между точками широты на экваторе больше, чем они вдоль полярных кругов, и это может быть решающий вопрос.

Моя цель состояла в том, чтобы разделить карту на сетку с шагом примерно 1 км по осям x и y. Таким образом, я взял минимальную и максимальную широту и долготу, вычислил их расстояния x и y от центра Стокгольма, а затем разделил диапазон широты и долготы на диапазон координат x и y (с использованием геосферы). Я сделал это, потому что хотел, чтобы точки были равномерно распределены при их построении (в противном случае расстояние между x точками вверху по сравнению с нижней частью карты было бы меньше из-за близости к экватору).

Затем я нанес эти точки на карту (с помощью ggmap) и заметил, что расстояние между точками в направлении y больше, чем в направлении x. Я полагаю, что карту можно было бы просто нарисовать в искаженном виде, но это кажется слишком искаженным, чтобы в это поверить. Я подозреваю, что делаю что-то не так, но не могу найти, что это может быть.

Пример кода ниже:

library("ggmap")
library("RgoogleMaps")
library("geosphere")

stockholm <- get_map("stockholm", zoom=11)
ggmap(stockholm)

places <- c('Tensta', 'Hanviken')

pos <- data.frame(Places = places, lat = NA, lon = NA, x = NA, y = NA)
reflatlon = getGeoCode('Stockholm, Sweden')

for(i in 1:length(places)) {
  latlon <- getGeoCode(paste0(places[i], ', Stockholm'))
  pos$lat[i] <- as.numeric(latlon[1])
  pos$lon[i] <- as.numeric(latlon[2])

  dist_y <- distGeo(c(latlon[1], reflatlon[2]), reflatlon) * sign(latlon[1] - reflatlon[1]) # same longitude
  dist_x <- distGeo(c(reflatlon[1], latlon[2]), reflatlon) * sign(latlon[2] - reflatlon[2]) # same latitude

  pos$x[i] <- dist_x
  pos$y[i] <- dist_y
}

deglatperm <- ( max(pos$lat) - min(pos$lat) ) / ( max(pos$y) - min(pos$y) ) # degrees latitude per metre
deglonperm <- ( max(pos$lon) - min(pos$lon) ) / ( max(pos$x) - min(pos$x) ) # degrees longitude per metre

seqlat <- seq(min(pos$lat), max(pos$lat), by = deglatperm*1000) # sequence with a point every ~1km
seqlon <- seq(min(pos$lon), max(pos$lon), by = deglonperm*1000) # sequence with a point every ~1km

seqlatlon <- expand.grid(seqlat, seqlon)
names(seqlatlon) <- c('lat', 'lon')

ggmap(stockholm) + geom_point(aes(x = lon, y = lat), data=seqlatlon)

выходной график

Как видно из выходного графика, расстояние между точками в направлении y по крайней мере в два раза больше, чем в направлении x.

Подводя итог: координаты x и y получены с использованием геосферы. Карта построена с использованием ggmap.

Я что-то не так делаю с геосферой? Или карты широты и долготы SO искажены? Когда я открываю Карты Google и использую инструмент «измерить расстояние» примерно между верхней и нижней, а также левой и правой точками, я получаю оценки в 16,3 и 16,9 км, тогда как значения, которые я получаю с геосферой, составляют 17 и 32 км (x и y ) соответственно.

Если бы кто-нибудь мог сказать мне, что здесь происходит, я был бы чрезвычайно благодарен!


person Granville Matheson    schedule 07.04.2016    source источник


Ответы (1)


Я стараюсь избегать работы в любой системе координат, где вместо расстояний используются градусы. В США я постоянно использую нашу систему State Plane. Судя по всему, Швеция использует систему RT. После того, как вы получите координаты вне градусов и на расстояния от точки отсчета, вы можете построить сетку, используя более обычные расстояния. Оттуда вы можете вернуть свои координаты в градусы, если хотите.

Я использую функцию spTranform для преобразования координат и использую функцию < руководство href = "http://spatialreference.org/ref/epsg/3022/" rel = "nofollow"> Пространственная справка, чтобы получить справочные коды для систем координат.

library("ggmap")
library("RgoogleMaps")
library("geosphere")
library("sp")

stockholm <- get_map("stockholm", zoom=11)
ggmap(stockholm)

places <- c('Tensta', 'Hanviken')

pos <- data.frame(Places = places, lat = NA, lon = NA)
reflatlon <- getGeoCode('Stockholm, Sweden')

for(i in 1:length(places)) {
  latlon <- getGeoCode(paste0(places[i], ', Stockholm'))
  pos$lat[i] <- as.numeric(latlon[1])
  pos$lon[i] <- as.numeric(latlon[2])
}

p <- SpatialPoints(data.frame(pos$lon, pos$lat), proj4string = CRS("+init=epsg:4326"))
p <- spTransform(p, CRS("+init=epsg:3022"))

seqx <- seq(min(p@coords[,1]), max(p@coords[,1]), by = 1000)
seqy <- seq(min(p@coords[,2]), max(p@coords[,2]), by = 1000)

pgrid <- expand.grid(seqx, seqy)
pgrid <- SpatialPoints(pgrid, proj4string = CRS("+init=epsg:3022"))
pgrid <- spTransform(pgrid, CRS("+init=epsg:4326"))
pgrid <- data.frame(pgrid@coords)

names(pgrid) <- c('lon', 'lat')

ggmap(stockholm) + geom_point(aes(x = lon, y = lat), data=pgrid)
person JMT2080AD    schedule 08.04.2016
comment
Отличный ответ - спасибо! Таким образом, проблема заключалась в том, что я использовал геосферу: я не осознавал, что системы координат такие сложные, но в ретроспективе это имеет смысл. Итак, если я правильно понимаю, вы преобразовали координаты в мировой системе координат epsg: 4362 в локальную топографическую систему координат epsg: 3022? И, выбирая именно epsg: 3022, вы сделали это, выполнив поиск на сайте spacereference.org? Я искал Стокгольм и не получил эту ссылку, но нашел 11 разных кодов. Или есть конкретный / автоматический способ / место для поиска того, какой код epsg использовать для указанного региона? Спасибо еще раз! - person Granville Matheson; 08.04.2016
comment
Преобразование в систему координат в единицах расстояния является ключевым. Чтобы найти систему координат RT90, я поискал в Google шведскую систему координат. Вы заметите, что на SR.org есть несколько систем координат RT. Сначала я попытался обыскать Стокгольм в СР, но тоже наткнулся на пустые руки. Я не знаю конкретного источника для поиска локальных систем координат (однако projfinder.org делает обратное, что довольно круто). EPSG - это просто номер по каталогу. Каждая система координат уникальна, и большинство из них признаны EPSG. Поиск в географических системах координат. Тяжелая штука. - person JMT2080AD; 08.04.2016
comment
Кроме того, вы можете отправить сообщение в GIS Stack Exchange, если ищете другие ответы. Если это работает для вас, вы можете проголосовать за него или отметить его как ответ. - person JMT2080AD; 08.04.2016
comment
Большой! Я искал в Google, чтобы попытаться понять это, но довольно быстро столкнулся с трудностями: действительно тяжелые вещи! Но я думаю, что это достаточно для меня, чтобы понять соответствующую систему координат в будущем, если я буду работать где-нибудь еще. Очень полезно! На самом деле это мой первый пост на StackOverflow, поэтому у меня недостаточно «репутации», чтобы проголосовать за него, но сейчас я помечу его как ответ - раньше не замечал эту кнопку. Спасибо еще раз! - person Granville Matheson; 08.04.2016
comment
Есть причина, по которой Интернет использует WGS 84 для всего - одна система координат, чтобы управлять ими всеми - я просто думаю, что это делает вычисления расстояний довольно сложными. Вероятно, есть инструмент, который будет правильно преобразовывать расстояние в градусы на лету. Опять же, обмен стеком ГИС может пролить больше света на эту проблему, если вам нужно сделать это для самых разных мест. Удачи! - person JMT2080AD; 08.04.2016