R：如何在定义 png 设备之前找出 OpenStreetMap 数据的尺寸（宽度，高度）？
R: How to find out dimension (width, height) of OpenStreetMap data before defining png device?
如何使用 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：
大圆距离在这里不起作用，因为图像是一块球面的变形表示（=矩形）。 （请原谅我的英语不好。）
大圆距离非常接近但不完全符合图像的纵横比。如果我从 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 传输的图像数据。
# 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 纵横比。
- 距离选择什么单位并不重要。大多数系统会让您选择英里、海里或公里。
所以我终于找到了解决方案。
我从 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.
相同
如何使用 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：
大圆距离在这里不起作用，因为图像是一块球面的变形表示（=矩形）。 （请原谅我的英语不好。）
大圆距离非常接近但不完全符合图像的纵横比。如果我从 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 传输的图像数据。
# 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 纵横比。
- 距离选择什么单位并不重要。大多数系统会让您选择英里、海里或公里。
所以我终于找到了解决方案。
我从 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.