Сглаживание карты ggplot2

Question

Предыдущие сообщения

Очистка карты с помощью geom_tile

Получить границы, чтобы пройти через штаты

Проблема / Вопрос

Я пытаюсь сгладить некоторые данные для сопоставления с ggplot2. Благодаря @MrFlick и @hrbrmstr я добился большого прогресса, но у меня возникают проблемы с получением эффекта "градиента" для перечисленных мне состояний.

Вот пример, чтобы дать вам представление о том, что я ищу:

**** Это именно то, чего я пытаюсь достичь.

http://nrelscience.org/2013/05/30/this-is-how-i-did-it-mapping-in-r-with-ggplot2/

(1) Как я могу максимально использовать ggplot2 с моими данными?

(2) Есть ли лучший метод для достижения эффекта градиента?

цели

Цели, которые я хотел бы достичь за счет этой награды:

(1) Интерполировать данные, чтобы построить растровый объект, а затем построить график с помощью ggplot2.

(или, если больше можно сделать с текущим графиком и растровый объект не является хорошей стратегией)

(2) Постройте лучшую карту с помощью ggplot2

Текущие результаты

Я играл со многими из этих разных графиков, но все еще не доволен результатами по двум причинам: (1) Градиент не говорит так много, как хотелось бы; и (2) презентация может быть улучшена, хотя я не уверен, как это сделать.

Как указывал @hrbrmstr, он мог бы дать лучшие результаты, если бы я выполнил некоторую интерполяцию с данными для получения большего количества данных, а затем поместил бы их в растровый объект и отобразил с помощью ggplot2. Я думаю, что это то, к чему я должен стремиться в данный момент, но я не уверен, как это сделать, учитывая мои данные.

Я перечислил ниже код, который я сделал до сих пор с результатами. Я действительно ценю любую помощь в этом вопросе. Спасибо.

Наборы данных

Вот два набора данных:

(1) Полный набор данных (175 МБ): PRISM_1895_db_all.csv (НЕ ДОСТУПНО)

https://www.dropbox.com/s/uglvwufcr6e9oo6/PRISM_1895_db_all.csv?dl=0

(2) Частичный набор данных (14 МБ): PRISM_1895_db.csv (НЕ ДОСТУПНО)

https://www.dropbox.com/s/0evuvrlm49ab9up/PRISM_1895_db.csv?dl=0

*** РЕДАКТИРОВАТЬ: Для тех, кто заинтересован, наборы данных не доступны, но я сделал пост на моем сайте, который связывает этот код с подмножеством данных Калифорнии на http://johnwoodill.com/pages/r-code.html

Участок 1

PRISM_1895_db <- read.csv("/.../PRISM_1895_db.csv")

regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas","minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")

ggplot() + 
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group)) +
  geom_point(data = PRISM_1895_db, aes(x = longitude, y = latitude, color = APPT), alpha = .5, size = 5) +
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group), color="white", fill=NA) +
  coord_equal()

Участок 2

PRISM_1895_db <- read.csv ("/.../ PRISM_1895_db.csv")

regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas","minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")

ggplot() + 
    geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group)) +
    geom_point(data = PRISM_1895_db, aes(x = longitude, y = latitude, color = APPT), alpha = .5, size = 5, shape = 15) +
    geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group), color="white", fill=NA) +
    coord_equal()

Участок 3

   PRISM_1895_db <- read.csv("/.../PRISM_1895_db.csv")

    regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas","minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")

ggplot() + 
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group)) +
  stat_summary2d(data=PRISM_1895_db, aes(x = longitude, y = latitude, z = APPT)) +
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group), color="white", fill=NA)

Answer 1

Пространственный вид CRAN заставил меня начать "Кригинг". Код ниже занимает ~7 минут для запуска на моем ноутбуке. Вы можете попробовать более простые интерполяции (например, некоторый вид сплайна). Вы также можете удалить некоторые местоположения из регионов с высокой плотностью. Вам не нужны все эти места, чтобы получить ту же карту тепла. Насколько я знаю, нет простого способа создать настоящий градиент с ggplot2 (gridSVG имеет несколько опций, но ничего похожего на "градиент сетки", который вы найдете в необычном редакторе SVG).

Как и требовалось, здесь есть интерполяция с использованием сплайнов (намного быстрее). Большая часть кода взята из черчения контуров на нерегулярной сетке.

Код для кригинга:

library(data.table)
library(ggplot2)
library(automap)

# Data munging
states=c("AR","IL","MO")
regions=c("arkansas","illinois","missouri")
PRISM_1895_db = as.data.frame(fread("./Downloads/PRISM_1895_db.csv"))
sub_data = PRISM_1895_db[PRISM_1895_db$state %in% states,c("latitude","longitude","APPT")]
coord_vars = c("latitude","longitude")
data_vars = setdiff(colnames(sub_data), coord_vars)
sp_points = SpatialPoints(sub_data[,coord_vars])
sp_df = SpatialPointsDataFrame(sp_points, sub_data[,data_vars,drop=FALSE])

# Create a fine grid
pixels_per_side = 200
bottom.left = apply(sp_points@coords,2,min)
top.right = apply(sp_points@coords,2,max)
margin = abs((top.right-bottom.left))/10
bottom.left = bottom.left-margin
top.right = top.right+margin
pixel.size = abs(top.right-bottom.left)/pixels_per_side
g = GridTopology(cellcentre.offset=bottom.left,
             cellsize=pixel.size,
             cells.dim=c(pixels_per_side,pixels_per_side))

# Clip the grid to the state regions
map_base_data = subset(map_data("state"), region %in% regions)
colnames(map_base_data)[match(c("long","lat"),colnames(map_base_data))] = c("longitude","latitude")
foo = function(x) {
  state = unique(x$region)
  print(state)
  Polygons(list(Polygon(x[,c("latitude","longitude")])),ID=state)
}
state_pg = SpatialPolygons(dlply(map_base_data, .(region), foo))
grid_points = SpatialPoints(g)
in_points = !is.na(over(grid_points,state_pg))
fit_points = SpatialPoints(as.data.frame(grid_points)[in_points,])

# Do kriging
krig = autoKrige(APPT~1, sp_df, new_data=fit_points)
interp_data = as.data.frame(krig$krige_output)
colnames(interp_data) = c("latitude","longitude","APPT_pred","APPT_var","APPT_stdev")

# Set up map plot
map_base_aesthetics = aes(x=longitude, y=latitude, group=group)
map_base = geom_polygon(data=map_base_data, map_base_aesthetics)
borders = geom_polygon(data=map_base_data, map_base_aesthetics, color="black", fill=NA)

nbin=20
ggplot(data=interp_data, aes(x=longitude, y=latitude)) + 
  geom_tile(aes(fill=APPT_pred),color=NA) +
  stat_contour(aes(z=APPT_pred), bins=nbin, color="#999999") +
  scale_fill_gradient2(low="blue",mid="white",high="red", midpoint=mean(interp_data$APPT_pred)) +
  borders +
  coord_equal() +
  geom_point(data=sub_data,color="black",size=0.3)

Код для сплайн-интерполяции:

library(data.table)
library(ggplot2)
library(automap)
library(plyr)
library(akima)

# Data munging
sub_data = as.data.frame(fread("./Downloads/PRISM_1895_db_all.csv"))
coord_vars = c("latitude","longitude")
data_vars = setdiff(colnames(sub_data), coord_vars)
sp_points = SpatialPoints(sub_data[,coord_vars])
sp_df = SpatialPointsDataFrame(sp_points, sub_data[,data_vars,drop=FALSE])

# Clip the grid to the state regions
regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas",
            "minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")
map_base_data = subset(map_data("state"), region %in% regions)
colnames(map_base_data)[match(c("long","lat"),colnames(map_base_data))] = c("longitude","latitude")
foo = function(x) {
  state = unique(x$region)
  print(state)
  Polygons(list(Polygon(x[,c("latitude","longitude")])),ID=state)
}
state_pg = SpatialPolygons(dlply(map_base_data, .(region), foo))

# Set up map plot
map_base_aesthetics = aes(x=longitude, y=latitude, group=group)
map_base = geom_polygon(data=map_base_data, map_base_aesthetics)
borders = geom_polygon(data=map_base_data, map_base_aesthetics, color="black", fill=NA)

# Do spline interpolation with the akima package
fld = with(sub_data, interp(x = longitude, y = latitude, z = APPT, duplicate="median",
                            xo=seq(min(map_base_data$longitude), max(map_base_data$longitude), length = 100),
                            yo=seq(min(map_base_data$latitude), max(map_base_data$latitude), length = 100),
                            extrap=TRUE, linear=FALSE))
melt_x = rep(fld$x, times=length(fld$y))
melt_y = rep(fld$y, each=length(fld$x))
melt_z = as.vector(fld$z)
level_data = data.frame(longitude=melt_x, latitude=melt_y, APPT=melt_z)
interp_data = na.omit(level_data)
grid_points = SpatialPoints(interp_data[,2:1])
in_points = !is.na(over(grid_points,state_pg))
inside_points = interp_data[in_points, ]

ggplot(data=inside_points, aes(x=longitude, y=latitude)) + 
  geom_tile(aes(fill=APPT)) + 
  stat_contour(aes(z=APPT)) +
  coord_equal() + 
  scale_fill_gradient2(low="blue",mid="white",high="red", midpoint=mean(inside_points$APPT)) +
  borders

Answer 2

Предыдущий ответ был, вероятно, не оптимальным (или точным) для ваших нужд. Это что-то вроде хака:

gg <- ggplot() 
gg <- gg + geom_polygon(data=subset(map_data("state"), region %in% regions), 
                        aes(x=long, y=lat, group=group))
gg <- gg + geom_point(data=PRISM_1895_db, aes(x=longitude, y=latitude, color=APPT), 
                      size=5, alpha=1/15, shape=19)
gg <- gg + scale_color_gradient(low="#023858", high="#ece7f2")
gg <- gg + geom_polygon(data=subset(map_data("state"), region %in% regions), 
                        aes(x=long, y=lat, group=group), color="white", fill=NA)
gg <- gg + coord_equal()
gg

это требует изменения size в geom_point для больших участков, но вы получите лучший эффект градиента, чем stat_summary2d поведение и передает ту же информацию.

Другим вариантом будет интерполировать больше APPT значения между имеющимися у вас значениями долготы и широты, затем преобразуйте их в более плотный растровый объект и нанесите его с помощью geom_raster как в приведенном вами примере.