R：如何在定义 png 设备之前找出 OpenStreetMap 数据的维度（宽度、高度）？答案

【问题标题】：R: How to find out dimension (width, height) of OpenStreetMap data before defining png device?R：如何在定义 png 设备之前找出 OpenStreetMap 数据的维度（宽度、高度）？
【发布时间】：2019-10-02 17:38:06
【问题描述】：

如何使用 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？我会想到……像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 传输的图像数据。

【问题讨论】：

您可以从左上角和右下角计算纵横比。
所以宽高比应该是(78-76)/(40-38) = 2/2 = 1。 2000x2000 示例表明情况并非如此。
# Washington DC upperLeft

标签： r png openstreetmap dimension

【解决方案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 以获得高/宽纵横比。
选择什么单位来表示距离并不重要。大多数系统都会让您选择英里、海里或公里。

【讨论】：

我会在最初的问题中给出更长的答案，因为 cmets 的空间有限。大圆距离非常接近，但并不完全符合图像的纵横比。无边框png文件的纵横比为0,778278278。用大圆计算的 NW->NE 与 NW->SW 的比率为0.7691904。 SW->SE 与 NE->SE 的比率更大，因为如果您靠近赤道，经度之间的距离会更大。它是0.7911577。但是，我可以将两个值的平均值作为近似值。但是我仍然缺乏最佳质量的总宽度和高度。
哦，对不起：我没看到你指的不是矩形的边缘而是中心线。我将这种情况添加到上面的问题中（补充 3）。这里的比例再次与 png 的纵横比不同：png = 0,778278278 vs. geo centerlines = 0.7802934。

【解决方案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 相同。

【讨论】：