ggplot/mapping графств США - проблемы с формами визуализации в R
Итак, у меня есть фрейм данных в R, называемый obesity_map, который в основном дает мне данные о состоянии, округе и уровне ожирения в округе. Это выглядит примерно так:
obesity_map = data.frame(state, county, obesity_rate)
Я пытаюсь визуализировать это на карте, показывая различные показатели ожирения в каждом округе США по следующей схеме:
us.state.map <- map_data('state')
head(us.state.map)
states <- levels(as.factor(us.state.map$region))
df <- data.frame(region = states, value = runif(length(states), min=0, max=100),stringsAsFactors = FALSE)
map.data <- merge(us.state.map, df, by='region', all=T)
map.data <- map.data[order(map.data$order),]
head(map.data)
map.county <- map_data('county')
county.obesity <- data.frame(region = obesity_map$state, subregion = obesity_map$county, value = obesity_map$obesity_rate)
map.county <- merge(county.obesity, map.county, all=TRUE)
ggplot(map.county, aes(x = long, y = lat, group=group, fill=as.factor(value))) + geom_polygon(colour = "white", size = 0.1)
И это в основном создает изображение, которое выглядит так:
Как вы можете видеть, США делятся на странные формы, цвета - это не один постоянный цвет в разных градиентах, и вы не можете многое из этого сделать. Но то, что я действительно хочу, это что-то вроде этого ниже, но с заполнением каждого графства:
Я довольно новичок в этом, поэтому я буду признателен за любую помощь!
Редактировать:
Вот вывод dput:
dput(obesity_map)
структура (список (X = 1:3141, FIPS = c(1L, 3L, 5L, 7L, 9L, 11L, 13L, 15L, 17L, 19L, 21L, 23L, 25L, 27L, 29L, 31L, 33L, 35L, 37л, 39л, 41л, 43л, 45л, 47л, 49л, 51л, 53л, 55л, 57л, 59л, 61л, 63л, 65л, 67л, 69л, 71л, 73л, 75л, 77л, 79л, 81л, 83л, 85л, 87л, 89л, 91л, 93л, 95л, 97л, 99л, 101л, 103л, 105л, 107л, 109л, 111л, 113л, 115л, 117л, 119л, 121л, 123л, 125л, 127л, 129л, 131л, 133л, 13л, 16л, 20л, 50л, 60л, 68л, 70л, 90л, 100л, 110л, 122л, 130л, 150л, 164л, 170л, 180л, 185л, 201л, 220л, 232л, 240л, 261л, 270л, 280л, 282л, 290л, 1л, 3л, 5л, 7л, 9л, 11л, 12л, 13л, 15л, 17л, 19л, 21л, 23л, 25л, 27л, 1л, 3л, 5л, 7л, 9л, 11л, 13л, 15л, 17л, 19L, 21L, 23L, 25L, 27L, 29L, 31L, 33L, 35L, 37L, 39L, 41L,
Это огромное количество цифр, потому что это для каждого округа США, поэтому я сократил результаты и вставил первые пару строк.
По сути, фрейм данных выглядит так:
print(head(obesity_map))
X FIPS state_names county_names obesity
1 1 1 Alabama Autauga 24.5
2 2 3 Alabama Baldwin 23.6
3 3 5 Alabama Barbour 25.6
4 4 7 Alabama Bibb 0.0
5 5 9 Alabama Blount 24.2
6 6 11 Alabama Bullock 0.0
Я также попытался использовать ggcounty, следуя приведенному примеру, но получаю сообщение об ошибке. Я не совсем уверен, что я сделал неправильно:
library(ggcounty)
# breaks
obesity_map$obese <- cut(obesity_map$obesity,
breaks=c(0, 5, 10, 15, 20, 25, 30),
labels=c("1", "2", "3", "4",
"5", "6"),
include.lowest=TRUE)
# get the US counties map (lower 48)
us <- ggcounty.us()
# start the plot with our base map
gg <- us$g
# add a new geom with our population (choropleth)
gg <- gg + geom_map(data=obesity_map, map=us$map,
aes(map_id=FIPS, fill=obesity_map$obese),
color="white", size=0.125)
Но я всегда получаю сообщение об ошибке: "Ошибка: аргумент должен быть принудительно приведен к неотрицательному целому числу"
Любая идея? Еще раз спасибо за вашу помощь! Я очень ценю это.
5 ответов
Может быть, немного поздно для другого ответа, но все же стоит поделиться, я думаю.
Чтение и предварительная обработка данных аналогичны ответу jlhoward с некоторыми отличиями:
library(tmap) # package for plotting
library(readxl) # for reading Excel
library(maptools) # for unionSpatialPolygons
# download data
download.file("http://www.ers.usda.gov/datafiles/Food_Environment_Atlas/Data_Access_and_Documentation_Downloads/Current_Version/DataDownload.xls", destfile = "DataDownload.xls", mode="wb")
df <- read_excel("DataDownload.xls", sheet = "HEALTH")
# download shape (a little less detail than in the other scripts)
f <- tempfile()
download.file("http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_us_050_00_20m.zip", destfile = f)
unzip(f, exdir = ".")
US <- read_shape("gz_2010_us_050_00_20m.shp")
# leave out AK, HI, and PR (state FIPS: 02, 15, and 72)
US <- US[!(US$STATE %in% c("02","15","72")),]
# append data to shape
US$FIPS <- paste0(US$STATE, US$COUNTY)
US <- append_data(US, df, key.shp = "FIPS", key.data = "FIPS")
Когда к объекту фигуры прикреплены правильные данные, хороплет можно нарисовать одной строкой кода:
qtm(US, fill = "PCT_OBESE_ADULTS10")
Это может быть улучшено путем добавления государственных границ, лучшей проекции и заголовка:
# create shape object with state polygons
US_states <- unionSpatialPolygons(US, IDs=US$STATE)
tm_shape(US, projection="+init=epsg:2163") +
tm_polygons("PCT_OBESE_ADULTS10", border.col = "grey30", title="") +
tm_shape(US_states) +
tm_borders(lwd=2, col = "black", alpha = .5) +
tm_layout(title="2010 Adult Obesity by County, percent",
title.position = c("center", "top"),
legend.text.size=1)
Так что это аналогичный пример, но пытается приспособить формат вашего obesity_map набор данных. Он также использует соединение таблиц данных, которое намного быстрее, чем merge(...)особенно с большими наборами данных, как у вас.
library(ggplot2)
# this creates an example formatted as your obesity.map - you have this already...
set.seed(1) # for reproducible example
map.county <- map_data('county')
counties <- unique(map.county[,5:6])
obesity_map <- data.frame(state_names=counties$region,
county_names=counties$subregion,
obesity= runif(nrow(counties), min=0, max=100))
# you start here...
library(data.table) # use data table merge - it's *much* faster
map.county <- data.table(map_data('county'))
setkey(map.county,region,subregion)
obesity_map <- data.table(obesity_map)
setkey(obesity_map,state_names,county_names)
map.df <- map.county[obesity_map]
ggplot(map.df, aes(x=long, y=lat, group=group, fill=obesity)) +
geom_polygon()+coord_map()
Кроме того, если в вашем наборе данных есть коды FIPS, что, как мне кажется, я настоятельно рекомендую вам использовать шейп-файл TIGER/Line Бюро переписей США (в котором также есть эти коды) и объединить его. Это гораздо надежнее. Например, при извлечении фрейма данных obesity_map регионы и регионы пишутся с большой буквы, тогда как во встроенном наборе данных округов в R их нет, поэтому вам придется с этим справиться. Кроме того, файл TIGER является актуальным, а внутренний набор данных - нет.
Так что это довольно интересный вопрос. Оказывается, что фактические данные о ожирении находятся на веб-сайте USDA и могут быть загружены здесь в виде файла MSExcel. На веб-сайте Бюро переписей здесь также есть шейп-файл округов США. И файл Excel, и шейп-файл содержат информацию FIPS. В R это может быть сравнительно просто:
library(XLConnect) # for loadWorkbook(...) and readWorksheet(...)
library(rgdal) # for readOGR(...)
library(RcolorBrewer) # for brewer.pal(...)
library(data.table)
setwd(" < directory with all your files > ")
wb <- loadWorkbook("DataDownload.xls") # from the USDA website
df <- readWorksheet(wb,"HEALTH") # this sheet has the obesity data
US.counties <- readOGR(dsn=".",layer="gz_2010_us_050_00_5m")
#leave out AK, HI, and PR (state FIPS: 02, 15, and 72)
US.counties <- US.counties[!(US.counties$STATE %in% c("02","15","72")),]
county.data <- US.counties@data
county.data <- cbind(id=rownames(county.data),county.data)
county.data <- data.table(county.data)
county.data[,FIPS:=paste0(STATE,COUNTY)] # this is the state + county FIPS code
setkey(county.data,FIPS)
obesity.data <- data.table(df)
setkey(obesity.data,FIPS)
county.data[obesity.data,obesity:=PCT_OBESE_ADULTS10]
map.df <- data.table(fortify(US.counties))
setkey(map.df,id)
setkey(county.data,id)
map.df[county.data,obesity:=obesity]
ggplot(map.df, aes(x=long, y=lat, group=group, fill=obesity)) +
scale_fill_gradientn("",colours=brewer.pal(9,"YlOrRd"))+
geom_polygon()+coord_map()+
labs(title="2010 Adult Obesity by Country, percent",x="",y="")+
theme_bw()
произвести это:
Это то, что я могу получить, работая с переменной отображения. Переименование его в "регион".
library(ggplot2)
library(maps)
m.usa <- map_data("county")
m.usa$id <- m.usa$subregion
m.usa <- m.usa[ ,-5]
names(m.usa)[5] <- 'region'
df <- data.frame(region = unique(m.usa$region),
obesity = rnorm(length(unique(m.usa$region)), 50, 10),
stringsAsFactors = F)
head(df)
region obesity
1 autauga 44.54833
2 baldwin 68.61470
3 barbour 52.19718
4 bibb 50.88948
5 blount 42.73134
6 bullock 59.93515
ggplot(df, aes(map_id = region)) +
geom_map(aes(fill = obesity), map = m.usa) +
expand_limits(x = m.usa$long, y = m.usa$lat) +
coord_map()
Я думаю, что все, что вам нужно было сделать, это переупорядочить переменную map.county, как вы делали это ранее для переменной map.data.
....
map.county <- merge(county.obesity, map.county, all=TRUE)
## reorder the map before plotting
map.county <- map.county[order(map.data$county),]
## plot
ggplot(map.county, aes(x = long, y = lat, group=group, fill=as.factor(value))) + geom_polygon(colour = "white", size = 0.1)
Опираясь на ответ @jlhoward: код с data.table терпит неудачу для меня таинственным образом:
Error in `:=`(FIPS, paste0(STATE, COUNTY)) :
Check that is.data.table(DT) == TRUE. Otherwise, := and `:=`(...) are defined for use in j, once only and in particular ways. See help(":=").
Эта ошибка случалась со мной несколько раз, но только тогда, когда код был внутри функции, даже просто минимальной оболочки. Это работало нормально в сценарии. Хотя теперь я не могу воспроизвести ошибку, я адаптировал его / ее код с merge() вместо data.table для полноты:
library(rgdal) # for readOGR(...)
library(ggplot2) # for fortify() and plot()
library(RColorBrewer) # for brewer.pal(...)
US.counties <- readOGR(dsn=".",layer="gz_2010_us_050_00_5m")
#leave out AK, HI, and PR (state FIPS: 02, 15, and 72)
US.counties <- US.counties[!(US.counties$STATE %in% c("02","15","72")),]
county.data <- US.counties@data
county.data <- cbind(id=rownames(county.data),county.data)
county.data$FIPS <- paste0(county.data$STATE, county.data$COUNTY) # this is the state + county FIPS code
df <- data.frame(FIPS=county.data$FIPS,
PCT_OBESE_ADULTS10= runif(nrow(county.data), min=0, max=100))
# Merge county.data to obesity
county.data <- merge(county.data,
df,
by.x = "FIPS",
by.y = "FIPS")
map.df <- fortify(US.counties)
# Merge the map to county.data
map.df <- merge(map.df, county.data, by.x = "id", by.y = "id")
ggplot(map.df, aes(x=long, y=lat, group=group, fill=PCT_OBESE_ADULTS10)) +
scale_fill_gradientn("",colours=brewer.pal(9,"YlOrRd"))+
geom_polygon()+coord_map()+
labs(title="2010 Adult Obesity by Country, percent",x="",y="")+
theme_bw()
Я немного новичок в использовании TMAP и пространственных данных, но подумал, что я отправлю в качестве продолжения Martijn Tennekes. Воспользовавшись его советом, я столкнулся с ошибкой на второй карте (с государственными границами). При запуске этой строки кода:
US_state <- unionSpatialPolygons(US,US$STATE)
Я продолжал получать эту ошибку: "Ошибка в unionSpatialPolygons(US, US$STATE): не объект SpatialPolygons"
Чтобы исправить, мне пришлось использовать другую переменную и запустить ее как фрейм данных пространственного полигона:
US <- read_shape("gz_2010_us_050_00_20m.shp")
US2<-readShapeSpatial("gz_2010_us_050_00_20m.shp")
US <- US[!(US$STATE %in% c("02","15","72")),]
US$FIPS <- paste0(US$STATE, US$COUNTY)
US <- append_data(US, med_inc_df, key.shp = "FIPS", key.data = "GEOID")
#the difference is here:
US_states <- unionSpatialPolygons(US2, US2$STATE)
tm_shape(US, projection="+init=epsg:2163") +
tm_polygons("estimate", border.col = "grey30", title="") +
tm_shape(US_states) +
tm_borders(lwd=2, col = "black", alpha = .5) +
tm_layout(title="2016 Median Income by County",
title.position = c("center", "top"),
legend.text.size=1)