R: Как определить размер (ширину, высоту) данных OpenStreetMap перед определением устройства png?

Question

Как я могу использовать OpenStreetMap для получения картографических данных и записи их в файл png в наилучшем качестве? Другими словами: Как узнать ширину и высоту данных OSM до определения устройства png?

Следующий пример получает данные OSM из Вашингтона. Я предполагаю, что размер в два раза явно не так:

1. 2000 х 2000 точек
2. 200 х 400 точек

#!/usr/bin/env Rscript 
library(OpenStreetMap) 
# Washington DC 
upperLeft <- c(40.00,-78.00) 
lowerRight <- c(38.00,-76.00) 
# get OpenStreetMap map 
map <- openmap(upperLeft, lowerRight, minNumTiles=4) 

# How to find out the correct size other than try&error??? 

pngWidth <- 2000 # => white border left and right 
pngHeight <- 2000 # => picture not sharp 
# open PNG device 
png("washingtondc_2000x2000.png", width=pngWidth, height=pngHeight) 
# avoid useless border 
par(mai=c(0,0,0,0)) # margin area in inches 
par(mar=c(0,0,0,0)) # margin area in number of lines (rows) of text 
par(xaxs="i", yaxs="i") # x- and y-axis won't be extended 
# plot on PNG device 
plot(map) 
# close PNG device 
dev.off() 

pngWidth <- 200 # => white border top and bottom 
pngHeight <- 400 # => text unreadable small 
# open PNG device 
png("washingtondc_0200x0400.png", width=pngWidth, height=pngHeight) 
# avoid useless border 
par(mai=c(0,0,0,0)) # margin area in inches 
par(mar=c(0,0,0,0)) # margin area in number of lines (rows) of text 
par(xaxs="i", yaxs="i") # x- and y-axis won't be extended 
# plot on PNG device 
plot(map) 
# close PNG device 
dev.off() 

Первая попытка (2000 х 2000) приводит в белых границ для слева и справа. Также он выглядит не резким, потому что пиксели карты OSM растянуты до слишком высоких 2000 пикселей.

Вторая попытка (200 х 400) приводит к границам в верхней и нижней и нечитаемый текст из-за слишком малым разрешением.

---

Дополнение 1:

Если я отрезать всю белую кайму из файла PNG 2000x2000 Оставшийся является 1555x1998. Соотношение сторон поэтому 1555/1998 = 0,778278278.

Теперь я беру пристальный взгляд на map <- openmap(upperLeft, lowerRight, minNumTiles=4):

> summary(map) 
     Length Class Mode 
tiles 1  -none- list 
bbox 2  -none- list 

> summary(map$tiles) 
    Length Class Mode 
[1,] 5  osmtile list 

> map$tiles 
[[1]] 
$colorData 
    [1] "#F1EEE8" "#F1EEE8" "#F1EEE8" "#F1EEE8" "#F1EEE8" "#F1EEE8" "#F4D5A7" 
<snip> 
[99996] "#DCDCDC" "#E6D5B7" "#F6C378" "#E98587" 
[ reached getOption("max.print") -- omitted 69159 entries ] 

$bbox 
<snip> 

$projection 
<snip> 

$xres 
[1] 466 

$yres 
[1] 363 

attr(,"class") 
[1] "osmtile" 

Здесь $xres и $yres выглядят интересно. Я не знаю, что означают эти данные, но я ожидал бы, что $xres будет меньше $yres. Во всяком случае, я снова разделил: $yres/$xres = 363/466 = 0,778969957. Это почти соотношение сторон от обрезанного изображения.

Означает ли это что-нибудь? Как я могу напрямую получить доступ к $xres? Я бы подумал, с. подобный map$tiles$xres но это NULL.

Я все еще в темноте. И одного аспекта недостаточно. Я бы хотел также иметь общую ширину и высоту для лучшего качества.

---

Дополнение 2:

@glenn-randers-pehrson:

Большое расстояние круг здесь не работает, потому что изображение деформированное представление (= прямоугольник) куска поверхности сферы. (Пожалуйста, извините мой плохой английский.)

Большие расстояния на расстоянии очень близки, но не совсем соответствуют соотношению сторон изображения. Если я отрезаю всю белую границу от файла PNX 2000x2000, оставшееся будет прямоугольником с 1555x1998 px. Соотношение сторон поэтому 1555/1998 = 0,778278278.

NW-> NE к NW-> соотношение SW вычисленный с большими окружностями 0.7691904:

library(geosphere) 
upperLeft <- c(40.00,-78.00) # lat lon 
lowerRight <- c(38.00,-76.00) # lat lon 
NW <- c(upperLeft[2], upperLeft[1]) # lon lat 
SW <- c(upperLeft[2], lowerRight[1]) # lon lat 
NE <- c(lowerRight[2], upperLeft[1]) # lon lat 
SE <- c(lowerRight[2], lowerRight[1]) # lon lat 

dist_NW_NE <- distGeo(NW, NE) 
dist_NW_SW <- distGeo(NW, SW) 
dist_NW_NE/dist_NW_SW 

[1] 0.7691904 

SW-> SE для NE-> SE отношение, конечно, больше, потому что расстояния между меридианами крупнее если вы ближе к экватору. Это 0.7911577:

dist_SW_SE <- distGeo(SW, SE) 
dist_NE_SE <- distGeo(NE, SE) 
dist_SW_SE/dist_NE_SE 

[1] 0.7911577 

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

---

Дополнение 3:

же расчет, как указано выше, но не с краями, но с осевыми. Отношение еще отличается от соотношения в формате PNG в: 0.7802934

N <- c(-77, 40) 
S <- c(-77, 38) 
W <- c(-78, 39) 
E <- c(-76, 39) 
dist_W_E <- distGeo(W, E) 
dist_N_S <- distGeo(N, S) 
dist_W_E/dist_N_S 

[1] 0.7802934 

Я думаю, что решение должно быть меньше на гео вычислений, но еще где-то в анализе данных изображения tranfered от OSM.

Answer 1

Итак, я наконец нашел решение.

Я получаю некоторые графические данные из OpenStreetMap:

library(OpenStreetMap) 
upperLeft <- c(40.00,-78.00) 
lowerRight <- c(38.00,-76.00) 
map <- openmap(upperLeft, lowerRight, minNumTiles=4) 

Теперь вопрос, какой ширины и длины этот образ есть? Ответ:

pngWidth <- map$tiles[[1]]$yres[1] 
pngHeight <- map$tiles[[1]]$xres[1] 

Это 363 x 466 px. Таким образом, PNG является: информация

fileName <- paste(c("washingtondc_", width, "x", height, ".png"), collapse='') 
png(fileName, width=pngWidth, height=pngHeight) 
plot(map) 
dev.off() 

---

Бонус: raster() дает хороший обзор:

library(OpenStreetMap) 
library(raster) 
raster <- raster(map) 
raster 

class  : RasterStack 
dimensions : 466, 363, 169158, 3 (nrow, ncol, ncell, nlayers) 
resolution : 613.3305, 614.8422 (x, y) 
extent  : -8682920, -8460281, 4579426, 4865942 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
names  : layer.1, layer.2, layer.3 
min values :  52,  51,  51 
max values :  253,  253,  253 

$xres из map$tiles путать меня, потому что он показывает значение высоты. Но raster() называет то же значение nrow, что имеет больше смысла, по крайней мере, для меня.То же самое для map$tiles[[1]]$yres[1] против raster()ncol.

Answer 2

# Washington DC 
upperLeft <- c(40.00,-78.00) 
lowerRight <- c(38.00,-76.00) 

Эти координаты находятся в градусах широты и долготы, а не в пикселях или милях.

Вы можете использовать бесплатную онлайн-версию GPS Visualizer service или подобное программное обеспечение для расчета большого расстояния между двумя точками.

• Введите 39, -76 и 39, -78, чтобы получить distanceEastWest
• Enter 38, -77 и 40, -77, чтобы получить distanceNorthSouth
• Разделить distanceNorthSouth от distanceEastWest, чтобы получить высота/ширина соотношение сторон.
• Неважно, какие единицы вы выбрали для расстояния. Большинство систем предложит вам выбрать мили, морские мили или километры.