Есть несколько примеров и руководств по составлению карт с использованием 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.
Что значит «R останавливается»? Сообщение об ошибке? Компьютер взрывается? Вы выполнили 'require (maptools)', правильно? Вы даже сделали 'install.packages (maptools)' да? Вы читали R Spatial TaskView? Если нет, найдите его и прочитайте. – Spacedman
Дополнительная информация. Я пытаюсь сначала прочитать шейп-файлы TIGER, надеясь объединить этот пространственный набор данных с моими данными и в конечном итоге построить. У меня возникла проблема в самом начале при попытке прочитать файл формы. Ниже приведен код с выходом --- Код --- требуется (maptools) фигуры <-readShapeSpatial ("tl_2013_us_zcta510.л.с ") --Output-- Ошибка: не может выделить вектор размером 317 Kb – vinod