6

使用O'Reilly 在 R 中的数据混搭作为灵感，我正在尝试在此处找到的犹他州盐湖县 shapefile 上绘制一些地址。

我有数据框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

所以看起来我已经拥有了绘制所需的一切——但它不起作用。当我加载 shapefile 并使用

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

我只能看着一个空的 shapefile。此外，当我尝试对我的 eventData 对象运行 findPolys 时，它返回 NULL。

> findPolys(addressEvents,myShapeFile)
NULL

我怎样才能使这项工作？我能够毫无问题地完成 O'Reilly 教程，并且很难弄清楚我在哪里出错了。我不知道它是 shapefile、我的数据框还是其他。

这是我用来导入数据和 shapefile 的命令

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)
4

2 回答 2

5

您可能还想查看这些相关问题，尤其是在 Eduardo 的回复中：

于 2009-09-26T19:03:01.707 回答
4

似乎 PBSmapping 使用一些粗略的启发式方法来计算 .prj 文件的投影。（请参阅帮助（importShapefile））。我个人不了解 prj 文件中的所有内容，但使用这个网站 www.spatialreference.org 我认为您的地图匹配

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

每当我得到一个新的形状文件时，我都会在这个网站上找到它的投影系统，然后查找 proj4 字符串，在这种情况下是 "+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=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"))

然后使用转换为 latlong

library(rgdal)

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

和

情节（sftransformed，轴= T）

给出坐标轴上具有正确单位的图.

不确定 PBSmapping 是否理解 proj4 字符串？好像不老实。

于 2009-09-26T13:27:33.177 回答