R：如何在定义 png 设备之前找出 OpenStreetMap 数据的尺寸（宽度，高度）？

Question

如何使用 OpenStreetMap 获取地图数据并将其写入质量最好的 png 文件？换句话说：如何在定义 png 设备之前知道 OSM 数据的宽度和高度？

以下示例从华盛顿特区获取 OSM 数据。我猜尺寸的两倍显然是错误的：

2000 x 2000 像素
200 x 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 x 2000) 导致左右出现白色边框。它看起来也不清晰，因为 OSM 地图像素被拉伸到太高的 2000 像素。

_{（来源：ibin.co）}

由于分辨率太小，第二次尝试 (200 x 400) 导致顶部和底部出现边框以及文本不可读。

_{（来源：ibin.co）}

加法 1：

如果我从 2000x2000 png 文件中切掉所有白色边框，剩下的就是 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？我本以为 s.th。像 map$tiles$xres 但那是 NULL.

我还是一头雾水。仅纵横比是不够的。我还想要最佳质量的总宽度和高度。

加法 2：

@glenn-randers-pehrson:

大圆距离在这里不起作用，因为图像是一块球面的变形表示（=矩形）。（请原谅我的英语不好。）

大圆距离非常接近但不完全符合图像的纵横比。如果我从 2000x2000 png 文件中切掉所有白色边框，剩下的就是一个 1555x1998 像素的矩形。因此纵横比为 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

我认为解决方案应该较少地进行地理计算，而更多地分析从 OSM 传输的图像数据。

Answer 1

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

这些坐标以纬度和经度为单位，而不是像素或英里。

您可以使用免费在线GPS Visualizer service或类似软件来计算两点之间的大圆距离。

输入 39,-76 和 39,-78 得到东西方距离
输入 38,-77 和 40,-77 得到南北距离
将 distanceNorthSouth 除以 distanceEastWest 得到 height/width 纵横比。
距离选择什么单位并不重要。大多数系统会让您选择英里、海里或公里。

Answer 2

所以我终于找到了解决方案。

我从 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 像素。所以 png 是：

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

_{（来源：ibin.co）}

奖金信息：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

map$tiles 中的 $xres 让我很困惑，因为它显示的是高度值。但是 raster() 命名相同的值 nrow 这至少对我来说更有意义。 map$tiles[[1]]$yres[1] 与 raster() 的 ncol.

相同

R：如何在定义 png 设备之前找出 OpenStreetMap 数据的尺寸（宽度，高度）？

R: How to find out dimension (width, height) of OpenStreetMap data before defining png device?

r

png

openstreetmap

dimension