3

你好 stackoverflow 社区！

有人可以帮助我解决我面临的 R GIS 问题。我正在尝试将一个可识别的变量分配给我已地理编码的地址列表。地理编码是从 Google Maps API 获得的，因此我有纬度和经度信息（即 -155.6019 18.99883）。我想将此信息与特定的形状文件一起使用。我的困境是我拥有的形状文件没有使用相同的纬度和经度系统。我附上了代码，以便您可以查看形状文件中使用的坐标系（即 843662.6 2132942）。

我想请教的是如何匹配我的地址列表和这个形状文件之间的坐标，以便我可以使用“叠加”功能将两者匹配在一起。

感谢您的时间！

#function to download shapefile from web
dlshapefile <- function(shploc,shpfile) {
temp <- tempfile()
download.file(shploc, temp)
unzip(temp)
}

temp <- tempfile()

require(maptools)
dlshapefile(shploc="http://files.hawaii.gov/dbedt/op/gis/data/highdist_n83.shp.zip", temp)
P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
shapeFile <- readShapePoly("highdist_n83.shp", verbose=TRUE, proj4string=P4S.latlon)
4

1 回答 1

7

我会readOGR从rgdal包中使用，因为它保留了投影信息。

require(rgdal)

in.dir <- "your_directory_name"

sh <- readOGR(in.dir, layer = "highdist_n83")
sh@proj4string # note projection details
sh2 <- spTransform(sh, CRS("+proj=longlat +datum=WGS84"))
sh2@proj4string # and after transforming
plot(sh2)

这给了我：

> sh@proj4string # note projection details
CRS arguments:
 +proj=utm +zone=4 +datum=NAD83 +units=m +no_defs +ellps=GRS80
+towgs84=0,0,0

和：

> sh2@proj4string # and after transforming
CRS arguments:
 +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

看看这是否适用于您要添加的地址数据。编辑如果你想沿着这ggplot2条路线走，你可以做这样的事情：

# Plotting with ggplot - first transform shapefile to data frame
sh2.df <- fortify(sh2)

# Fake address data taken from Google maps
add <- data.frame(name = c("University of Hawaii", "Mokuleia Beach Park"),
                      lat = c(21.298971, 21.580508),
                      long = c(-157.817722, -158.191017))

# And plot
ggplot(data = sh2.df, aes(x = long, y = lat, group = group)) +
     geom_polygon(colour = "grey20", fill = "white", size = 0.8) +
     geom_point(data = add, size = 3, aes(x = long, y = lat, group = NULL), pch = 21, fill = "red") +
     coord_equal() +
     theme()

这给出了这个：

截屏

于 2013-06-01T05:32:05.463 回答