2013-12-25 6 views
4

Есть вопрос по картированию с R, в частности, вокруг фоновой картограммы в R.фоновая картограмма в R - ТИГР Shapefile вопроса

у меня есть набор данные почтовых индексов, присвоенные есть и некоторые связанные с ним данные (набор данные here).

Моим окончательным форматом данных является: ID области, ZIP, Вероятностная ценность, количество клиентов, вероятность области и общая сумма клиентов. Я пытаюсь представить эти данные, построив вероятность области и Area Customer Total на карте. Я попытался сделать это, используя TIGER Shapefiles, но я думаю, что R не может справиться с полной страной.

Мне комфортно со статистическими возможностями, и теперь я перемещаю все свои карты с сторонних приложений, ориентированных на ГИС, на выполнение всех моих карт в R. Есть ли у кого-нибудь какие-либо указания на то, как достичь этого изнутри?

Чтобы быть немного более подробно, вот та точка, где R перестает работать -

shapes <- readShapeSpatial("tl_2013_us_zcta510.shp") 

(где файл л.с является/Перепись TIGER) файл форма.

Редактировать - Предоставляются дальнейшие детали. Я пытаюсь сначала прочитать шейп-файлы TIGER, надеясь объединить этот пространственный набор данных с моими данными и в конечном итоге построить. У меня возникла проблема в самом начале при попытке прочитать файл формы. Ниже приведен код, с выходом

require(maptools) 
shapes<-readShapeSpatial("tl_2013_us_zcta510.shp") 

Error: cannot allocate vector of size 317 Kb 
+0

Что значит «R останавливается»? Сообщение об ошибке? Компьютер взрывается? Вы выполнили 'require (maptools)', правильно? Вы даже сделали 'install.packages (maptools)' да? Вы читали R Spatial TaskView? Если нет, найдите его и прочитайте. – Spacedman

+0

Дополнительная информация. Я пытаюсь сначала прочитать шейп-файлы TIGER, надеясь объединить этот пространственный набор данных с моими данными и в конечном итоге построить. У меня возникла проблема в самом начале при попытке прочитать файл формы. Ниже приведен код с выходом --- Код --- требуется (maptools) фигуры <-readShapeSpatial ("tl_2013_us_zcta510.л.с ") --Output-- Ошибка: не может выделить вектор размером 317 Kb – vinod

ответ

7

Есть несколько примеров и руководств по составлению карт с использованием R, но большинство из них очень общего характера и, к сожалению, большинство картографических проектов есть нюансы, которые создают непостижимые проблемы. Это ваш пример.

Самая большая проблема, с которой я столкнулся, заключалась в том, что форма шейпинга области индекса почтового индекса США для всего США огромна: ~ 800 МБ. При загрузке с использованием readOGR(...) объект R SpatialPolygonDataFrame составляет около 913 МБ. Попытка обработать файл такого размера (например, преобразование в фрейм данных с использованием fortify(...)), по крайней мере в моей системе, привело к ошибкам, подобным указанным выше. Таким образом, решение заключается в подмножестве файла, основанного на почтовых индексах, которые фактически находятся в ваших данных.

Эта карта:

был сделан на основе данных, используя следующий код.

library(rgdal) 
library(ggplot2) 
library(stringr) 
library(RColorBrewer) 

setwd("<directory containing shapfiles and sample data>") 

data  <- read.csv("Sample.csv",header=T) # your sample data, downloaded as csv 
data$ZIP <- str_pad(data$ZIP,5,"left","0") # convert ZIP to char(5) w/leading zeros 

zips  <- readOGR(dsn=".","tl_2013_us_zcta510") # import zip code polygon shapefile 
map  <- zips[zips$ZCTA5CE10 %in% data$ZIP,] # extract only zips in your Sample.csv 
map.df <- fortify(map)  # convert to data frame suitable for plotting 
# merge data from Samples.csv into map data frame 
map.data <- data.frame(id=rownames([email protected]),[email protected]$ZCTA5CE10) 
map.data <- merge(map.data,data,by="ZIP") 
map.df <- merge(map.df,map.data,by="id") 
# load state boundaries 
states <- readOGR(dsn=".","gz_2010_us_040_00_5m") 
states <- states[states$NAME %in% c("New York","New Jersey"),] # extract NY and NJ 
states.df <- fortify(states) # convert to data frame suitable for plotting 

ggMap <- ggplot(data = map.df, aes(long, lat, group = group)) 
ggMap <- ggMap + geom_polygon(aes(fill = Probability_1)) 
ggMap <- ggMap + geom_path(data=states.df, aes(x=long,y=lat,group=group)) 
ggMap <- ggMap + scale_fill_gradientn(name="Probability",colours=brewer.pal(9,"Reds")) 
ggMap <- ggMap + coord_equal() 
ggMap 

Объяснение:

Пакет rgdal облегчает создание R пространственных объектов из ESRI шейпфайлов. В вашем случае мы импортируем файл формы polygon в объект SpatialPolygonDataFrame в R. Последний имеет две основные части: секцию многоугольника, которая содержит точки широты и долготы, которые будут соединены для создания многоугольников на карте, и раздел данных который содержит информацию о многоугольниках (так, одна строка для каждого многоугольника). Если, например, мы вызываем объект Spatial map, то на два раздела можно ссылаться как [email protected] и [email protected]. Основная проблема при создании choropleth карт - связать данные из вашего файла Sample.csv с соответствующими полигонами (почтовые индексы).

Так основной рабочий процесс выглядит следующим образом:

1. Load polygon shapefiles into Spatial object (=> zips) 
2. Subset if appropriate (=> map). 
3. Convert to data frame suitable for plotting (=> map.df). 
4. Merge data from Sample.csv into map.df. 
5. Draw the map. 

Шаг 4 является тот, который вызывает все проблемы. Сначала мы должны сопоставлять почтовые коды с каждым полигоном. Затем мы должны связать Probability_1 с каждым почтовым индексом. Это трехэтапный процесс.

Каждый многоугольник в файле пространственных данных имеет уникальный идентификатор, но эти идентификаторы не являются почтовыми индексами. Идентификаторы многоугольников хранятся в именах строк в [email protected]. Почтовые индексы хранятся в [email protected], в столбце ZCTA5CE10. Поэтому сначала мы должны создать кадр данных, который связывает имена строк [email protected] (id) с [email protected]$ZCTA5CE10 (ZIP). Затем мы объединяем ваш Sample.csv с результатом, используя поле ZIP в обоих кадрах данных. Затем мы объединим результат этого в map.df. Это можно сделать в 3 строках кода.

Рисование карты включает в себя рассказ ggplot о том, какой набор данных использовать (map.df), какие столбцы использовать для x и y (long и lat) и как сгруппировать данные по многоугольнику (group = group). Колонки long, lat и group в map.df созданы по вызову fortify(...). Вызов geom_polygon(...) указывает ggplot на рисование полигонов и заполнение с использованием информации в map.df$Probability_1. Вызов geom_path(...) сообщает ggplot о создании слоя с границами состояний. Вызов scale_fill_gradientn(...) сообщает ggplot использовать цветовую схему, основанную на палитре цветного пивовара «Красные». Наконец, вызов coord_equal(...) сообщает ggplot использовать тот же масштаб для x и y, чтобы карта не искажалась.

NB: Государственный пограничный слой использует US States TIGER file.

+0

Спасибо Это полезно – vinod

2

Я бы посоветовал следующее.

  • Использование readOGR из rgdal пакета, а не readShapeSpatial.
  • Рассмотрите возможность использования ggplot2 для красивейших карт - многие из примеров используют это.
  • Обратитесь к одному из существующих примеров создания choropleth, такого как this one, чтобы получить обзор.
  • Начните с простого choropleth и постепенно добавьте свои данные; не пытайтесь сразу все исправить.
  • Если вам нужна дополнительная справка, создайте reproducible example с помощью SMALL fake dataset и со ссылками на соответствующие шейп-файлы. Идея заключается в том, что вы упрощаете нам помощь, а не препятствуете нам, не поставляя код и данные в ваш вопрос.
+0

Ответ ниже гораздо больше заслуживает того, чтобы быть выбран в качестве правильного ответа, чем мой (это один).. Если вы согласитесь, вы должны уметь отменить этот ответ, нажав на отметку, оставив вас свободным, чтобы выбрать другой ответ в качестве правильного. – SlowLearner

+0

Спасибо, это прекрасно. – vinod