Помогите построить географические данные в R с помощью PBSMapping и шейп-файлов

Используя в качестве вдохновения Мэшапы данных О'Рейли в 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, он возвращает NULL.

> findPolys(addressEvents,myShapeFile)
NULL

Как я могу заставить это работать? Я смог пройти обучение O'Reilly без каких-либо проблем, и мне сложно понять, в чем я ошибаюсь. Я не знаю, это шейп-файл, мой фрейм данных или что-то еще.

Вот команды, которые я использую для импорта данных и шейп-файла.

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)

person Neil Kodner    schedule 26.09.2009    source источник
comment
Мое чутье подсказывает мне, что это не мои данные, а мой шейп-файл. Я новичок в концепции шейп-файлов. Когда я рисую Polys, используя пример шейп-файла О'Рейли, оси X и Y выглядят как длинная и широта. Когда я рисую Polys с использованием шейп-файлов штата Юта, оси X и Y выглядят как другая система нумерации.   -  person Neil Kodner    schedule 26.09.2009
comment
Я никогда не пользовался пакетом PBSmapping. Я немного испортил maptools.   -  person JD Long    schedule 28.09.2009


Ответы (2)


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

person Shane    schedule 26.09.2009

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

и

сюжет (преобразованный, оси = T)

дает сюжет с правильными узлами по осям.

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

person Ira Cooke    schedule 26.09.2009