帮助绘制地理数据R使用PBSMapping和形状文件

Question

使用 O'reilly的数据混合在R 作为灵感，我想对积极少数的地址在一个要素的盐湖县，犹他州发现 在这里，.

我有数据的框架geoTable:

> geoTable
         address        Y         X EID
1    130 E 300 S 40.76271 -111.8872   1
2    875 E 900 S 40.74992 -111.8660   2
3   2200 S 700 E 40.72298 -111.8714   3
4    702 E 100 S 40.76705 -111.8707   4
5 177 East 200 S 40.76518 -111.8859   5
6    702 3rd ave 40.77264 -111.8683   6
7   2175 S 900 E 40.72372 -111.8652   7
8   803 E 2100 S 40.72556 -111.8680   8

我已经迫使它变成一个eventData对象：

> addressEvents<-as.EventData(geoTable,projection=NA)
> addressEvents
         address        Y         X EID
1    130 E 300 S 40.76271 -111.8872   1
2    875 E 900 S 40.74992 -111.8660   2
3   2200 S 700 E 40.72298 -111.8714   3
4    702 E 100 S 40.76705 -111.8707   4
5 177 East 200 S 40.76518 -111.8859   5
6    702 3rd ave 40.77264 -111.8683   6
7   2175 S 900 E 40.72372 -111.8652   7
8   803 E 2100 S 40.72556 -111.8680   8

因此它看起来像我们有我需要的一切情节，但它没有工作。当我加载的要素和情节使用

addPoints(addressEvents,col="red",cex=.5)

我离开寻找一个空的要素.此外，当我尝试和运行findPolys对我eventData目的，它返回空。

> findPolys(addressEvents,myShapeFile)
NULL

我如何可以使这项工作？我能够完成的O'reilly的教程中没有任何问题，我有困难搞清楚我要去哪里错在这里。我不知道，如果它的要素，我的数据框架，或whateverelse.

这里是命令我使用我的数据和形状

slc<-read.table('~/utah.txt',sep=',',header=TRUE,strip.white=TRUE,stringsAsFactors=FALSE)

myShapeFile<-importShapefile("/Users/neil/Downloads/SGID93_DEMOGRAPHIC_CensusTracts2000/SGID93_DEMOGRAPHIC_CensusTracts2000",readDBF=TRUE)

Answer 1

你也可能想看看这些相关的问题，尤其是在爱德华多*的答复:

• 绘色地图-with-zip-码-在-r-或蟒蛇
• 发展中地理的专题地图与-r

Answer 2

好像PBSmapping使用一些粗试探法来从.prj文件制定出的投影。 （见帮助（importShapefile））。我个人不明白PRJ文件里面的所有东西，但在使用本网站www.spatialreference.org我看你的地图匹配

http://www.spatialreference.org/ref/epsg/26912/

每当我得到一个新的形状文件，我觉得这是投影系统本网站上，然后查找proj4字符串，在这种情况下是 “+ PROJ = UTM +区= 12 + ellps = GRS80 +基准= NAD83 +单位= M + no_defs”

（就像我说的，我不知道PBSmapping，但你可以使用maptools如下阅读本）

library(maptools)
sf=readShapeSpatial("SGID93_DEMOGRAPHIC_CensusTracts2000.shp",proj4string=CRS("+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

，然后将转换到使用经纬度

library(rgdal)

sftransformed=spTransform(sf,CRS("+proj=longlat"))

和

积（sftransformed，轴= T）

给出了在轴的右边单元的曲线图。

不知道PBSmapping理解一个proj4字符串？看起来不诚实。