Ajude a plotar dados geográficos em r usando PBSMapping e Shapefiles

Question

Usando Os mashups de dados de O'Reilly em r Como inspiração, estou tentando traçar um punhado de endereços em um shapefile de Salt Lake County, Utah encontrado aqui.

Eu tenho o quadro de dados geotável:

> 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

E eu o coagi a um objeto 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

Parece que tenho tudo o que preciso para plotar, mas não está funcionando. Quando carrego o shapfile e o enredo usando

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

Fiquei olhando para um shapefile vazio. Além disso, quando tento executar o FindPolys no meu objeto EventData, ele retorna nulo.

> findPolys(addressEvents,myShapeFile)
NULL

Como posso fazer isso funcionar? Consegui concluir o tutorial de O'Reilly sem problemas e estou tendo dificuldades para descobrir onde estou dando errado aqui. Não sei se é o shapefile, meu quadro de dados ou qualquer outra coisa.

Aqui estão os comandos que eu uso para importar meus dados e 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)

Answer 1

Você também pode querer olhar para essas perguntas relacionadas, especialmente nas respostas de Eduardo:

• plotar-cor-map-with-zip-codos-in-r-or-python
• Mapas-geográficos-geográficos em desenvolvimento

Answer 2

Parece que o PBSMapping usa algumas heurísticas brutas para elaborar a projeção do arquivo .prj. (Veja Ajuda (ImportshapeFile)). Pessoalmente, não entendo todas as coisas dentro de um arquivo PRJ, mas usando este site www.spatialreference.org eu acho que seu mapa corresponde

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

Sempre que recebo um novo arquivo de forma, acho seu sistema de projeção neste site e procure a string proj4, que neste caso é " +proj = utm +zone = 12 +ellps = grs80 +datum = nad83 +unidades = m + no_defs "

(Como eu disse, não conheço PBSMapping, mas você pode ler isso usando o Maptools da seguinte maneira)

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

e depois converta em latlongo usando

library(rgdal)

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

e

plot (sftransformed, eixos = t)

dá um gráfico com as unidades corretas nos eixos.

Não tenho certeza se o PBSMapping entende uma string proj4? Parece que não é para ser honesto.