在 shapefile 中保存一个 ggplot2 coord_map() 图表
Save a ggplot2 coord_map() chart in shapefile
我需要在 CARTO(又名 cartodb)中导出等高线图,因此我试图将此 stat2density 图表保存为地理数据文件格式,如 shapefile 或 geojson。
我可以使用 ggsave 将其保存为 SVG,但将其转换为 spdf 或 sf oblejct 将非常有帮助。
library(ggplot2)
library(ggmap)
data("crime")
crime<- head(crime,1000)
gg <- ggplot(aes(x = lon, y = lat), data=crime) +
stat_density2d(aes(alpha=..level.., color=..level.., fill=..level..),geom='polygon', bins = 10, size=0.5) +
scale_color_gradient(low = "grey", high = "#444444", guide = F)+
scale_fill_gradient(low = "yellow", high = "red", guide = F)+
scale_alpha( guide = F)+
coord_map()+
ggthemes::theme_map()
有什么想法吗?
这是一个解决方案,它结合了上述@hrbrmstr 和@Rich Pauloo 提出的想法,以及:
步骤 1。从ggplot对象中提取相关数据:
library(dplyr)
# return a list of data frames (each data frame contains coordinates for one contour);
# note that there may be multiple contours at the same alpha / colour / fill,
# hence the need to split by group rather than by these aesthetic mappings.
dg <- layer_data(gg) %>%
select(group, x, y) %>%
split(.$group) %>%
lapply(function(d){d[,-1]})
步骤 2。将数据框转换为 SpatialPolygonsDataFrame 对象,传递给 writeOGR
:
library(sp)
# convert each data frame to a Polygon class object
polygons <- lapply(dg, Polygon)
# convert each Polygon class object to Polygons class object
polygons <- lapply(seq_along(polygons),
function(i){
Polygons(list(polygons[[i]]),
ID = names(dg)[i])
})
# convert list of Polygons class object to one SpatialPolygons object
polygons <- SpatialPolygons(polygons)
# convert SpatialPolygons object to SpatialPolygonsDataFrame object
polygons <- SpatialPolygonsDataFrame(polygons,
data = layer_data(gg) %>%
select(group, alpha, colour, fill) %>%
unique(),
match.ID = "group")
步骤 3。将 SpatialPolygonsDataFrame 对象保存为 shapefile:
rgdal::writeOGR(obj = polygons,
dsn = getwd(), # or wherever you wish to store it
layer = "filename", # or whatever you wish to name it
driver = "ESRI Shapefile")
在 R 中验证结果(我更愿意在单独的 GIS 程序中验证它,但我没有在这台计算机上安装任何程序):
# read the shapefile back into R
sp <- rgdal::readOGR(dsn = getwd(),
layer = "filename")
# fortify as a data frame
spdf <- left_join(broom::tidy(sp, region = "group"),
sp@data,
by = c("id" = "group"))
# plot
ggplot(spdf,
aes(x = long, y = lat, group = group, alpha = alpha)) +
geom_polygon(color = "black") +
coord_map()
我需要在 CARTO(又名 cartodb)中导出等高线图,因此我试图将此 stat2density 图表保存为地理数据文件格式,如 shapefile 或 geojson。 我可以使用 ggsave 将其保存为 SVG,但将其转换为 spdf 或 sf oblejct 将非常有帮助。
library(ggplot2)
library(ggmap)
data("crime")
crime<- head(crime,1000)
gg <- ggplot(aes(x = lon, y = lat), data=crime) +
stat_density2d(aes(alpha=..level.., color=..level.., fill=..level..),geom='polygon', bins = 10, size=0.5) +
scale_color_gradient(low = "grey", high = "#444444", guide = F)+
scale_fill_gradient(low = "yellow", high = "red", guide = F)+
scale_alpha( guide = F)+
coord_map()+
ggthemes::theme_map()
有什么想法吗?
这是一个解决方案,它结合了上述@hrbrmstr 和@Rich Pauloo 提出的想法,以及
步骤 1。从ggplot对象中提取相关数据:
library(dplyr)
# return a list of data frames (each data frame contains coordinates for one contour);
# note that there may be multiple contours at the same alpha / colour / fill,
# hence the need to split by group rather than by these aesthetic mappings.
dg <- layer_data(gg) %>%
select(group, x, y) %>%
split(.$group) %>%
lapply(function(d){d[,-1]})
步骤 2。将数据框转换为 SpatialPolygonsDataFrame 对象,传递给 writeOGR
:
library(sp)
# convert each data frame to a Polygon class object
polygons <- lapply(dg, Polygon)
# convert each Polygon class object to Polygons class object
polygons <- lapply(seq_along(polygons),
function(i){
Polygons(list(polygons[[i]]),
ID = names(dg)[i])
})
# convert list of Polygons class object to one SpatialPolygons object
polygons <- SpatialPolygons(polygons)
# convert SpatialPolygons object to SpatialPolygonsDataFrame object
polygons <- SpatialPolygonsDataFrame(polygons,
data = layer_data(gg) %>%
select(group, alpha, colour, fill) %>%
unique(),
match.ID = "group")
步骤 3。将 SpatialPolygonsDataFrame 对象保存为 shapefile:
rgdal::writeOGR(obj = polygons,
dsn = getwd(), # or wherever you wish to store it
layer = "filename", # or whatever you wish to name it
driver = "ESRI Shapefile")
在 R 中验证结果(我更愿意在单独的 GIS 程序中验证它,但我没有在这台计算机上安装任何程序):
# read the shapefile back into R
sp <- rgdal::readOGR(dsn = getwd(),
layer = "filename")
# fortify as a data frame
spdf <- left_join(broom::tidy(sp, region = "group"),
sp@data,
by = c("id" = "group"))
# plot
ggplot(spdf,
aes(x = long, y = lat, group = group, alpha = alpha)) +
geom_polygon(color = "black") +
coord_map()