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

Question

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

以下示例从华盛顿特区获取OSM数据。我猜测大小明显是错误的两次：

1. 2000 x 2000 px
2. 200 x 400 px

#!/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 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

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

Answer 1

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

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

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

• 输入39，-76和39，-78来获取EastEast的距离
• 输入38，-77和40，-77获得北北距离
• 将distanceNorthSouth除以distanceEastWest以得到高度/宽度长宽比。
• 选择距离的单位并不重要。大多数系统将为您提供英里，海里或公里的选择。

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 

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