使用 R 在 rgdal 包中的地图上绘制点
Plotting points on a map in the rgdal package with R
借用 Rob Berry (http://rob-barry.com/2015/06/14/Mapping-in-R/) 的代码,我制作了一张纽约市地图。我想在地图上绘制许多经纬度点。问题是绘制这样的地图会导致绘图区域超出合理的经纬度范围,所以我认为必须有一种方法可以将我的点转换为地图比例,或者重新缩放地图,以便绘图 space 可以使用 points() 函数放置经纬度点。
这是 Rob Berry 的代码:
download.file(destfile = "nypp_15b.zip")
unzip(zipfile = "nypp_15b.zip")
library("rgdal")
nypp <- readOGR("nypp_15b", "nypp")
plot(nypp)
漂亮的地图!但现在请注意绘图范围:
par(“usr”)
地块 space 数字看起来像 888196.7、1092361.0、114013.0、278953.2,所以很明显像下面这些经纬度点不会显示在地图上。那么如何让我的点在地图上正确绘制?
lat <- c(40.75002, 40.74317)
lon <- c(-73.96905 -74.00366)
以下内容不起作用,因为比例如此不同:
points(lat,lon, col = “red”)
非常感谢。
nypp 在投影坐标系中,因此您需要更改为您的点或 nypp 的坐标系。你可以这样做:-
nypp <- readOGR("nypp_15b", "nypp")
## Check the CRS of nypp
crs(nypp)
## CRS arguments:
+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000
+y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
plot(nypp)
lat <- c(40.75002, 40.74317)
lon <- c(-73.96905, -74.00366)
df <- data.frame(lat, lon)
## Convert to spatial dataframe and change the coordinates of the points
coordinates(df) <- ~lon + lat
crs(df) <- CRS("+proj=longlat +datum=WGS84")
df <- spTransform(df, crs(nypp))
## Add points to the plot
points(df$lon, df$lat, col = "red", pch =19)
结果:
借用 Rob Berry (http://rob-barry.com/2015/06/14/Mapping-in-R/) 的代码,我制作了一张纽约市地图。我想在地图上绘制许多经纬度点。问题是绘制这样的地图会导致绘图区域超出合理的经纬度范围,所以我认为必须有一种方法可以将我的点转换为地图比例,或者重新缩放地图,以便绘图 space 可以使用 points() 函数放置经纬度点。 这是 Rob Berry 的代码:
download.file(destfile = "nypp_15b.zip")
unzip(zipfile = "nypp_15b.zip")
library("rgdal")
nypp <- readOGR("nypp_15b", "nypp")
plot(nypp)
漂亮的地图!但现在请注意绘图范围:
par(“usr”)
地块 space 数字看起来像 888196.7、1092361.0、114013.0、278953.2,所以很明显像下面这些经纬度点不会显示在地图上。那么如何让我的点在地图上正确绘制?
lat <- c(40.75002, 40.74317)
lon <- c(-73.96905 -74.00366)
以下内容不起作用,因为比例如此不同:
points(lat,lon, col = “red”)
非常感谢。
nypp 在投影坐标系中,因此您需要更改为您的点或 nypp 的坐标系。你可以这样做:-
nypp <- readOGR("nypp_15b", "nypp")
## Check the CRS of nypp
crs(nypp)
## CRS arguments:
+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000
+y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
plot(nypp)
lat <- c(40.75002, 40.74317)
lon <- c(-73.96905, -74.00366)
df <- data.frame(lat, lon)
## Convert to spatial dataframe and change the coordinates of the points
coordinates(df) <- ~lon + lat
crs(df) <- CRS("+proj=longlat +datum=WGS84")
df <- spTransform(df, crs(nypp))
## Add points to the plot
points(df$lon, df$lat, col = "red", pch =19)
结果: