Помогите построить географические данные в R с использованием PBSMapping и Shapefiles

Question

С использованием О'Рейли Как вдохновение, я пытаюсь построить несколько адресов на форме -файле округа Солт -Лейк, штат Юта, найдена здесь.

У меня есть кадр данных 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

Похоже, что у меня есть все, что мне нужно, но это не работает. Когда я загружаю форму Shape и график, используя

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

Я остался смотреть на пустой форму для формы. Кроме того, когда я пытаюсь запустить FindPolys против моего объекта EventData, он возвращает NULL.

> findPolys(addressEvents,myShapeFile)
NULL

Как я могу заставить это работать? Я смог выполнить учебник O'Reilly без каких -либо проблем, и мне трудно понять, где я ошибаюсь здесь. Я не знаю, является ли это ShapeFile, My Data Frame или What Elelse.

Вот команды, которые я использую для импорта моих данных и 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

Вы также можете посмотреть на эти связанные вопросы, особенно в ответах Эдуардо:

• Построение цветов-карты-zip-codes-in-r-or-python
• развивающиеся географические тематические карты-с-R

Answer 2

Похоже, что 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"))

а затем конвертируйте в Латлонг, используя

library(rgdal)

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

а также

График (sftransformed, Axes = t)

дает график с правыми единицами на оси.

Не уверен, понимает ли PBSMapping строку ProJ4? Похоже, это не честно.