R语言整理城市级碳排放数据并绘制出图 |
您所在的位置：网站首页 › 中国地图坐标线 › R语言整理城市级碳排放数据并绘制出图 |
R语言整理城市级碳排放数据并绘制出图
|
双碳这几年很火，笔者找到一份城市级别比较详细的碳排放数据，来自中国城市温室气体工作组，具体链接见： 城市温室气体 (cityghg.com)，数据以excel格式对公众开放，原始数据没有空间信息（经纬度坐标等），对于一个具有强迫症的GISer来说简直要命，就趁今天（2023.3.22）有空搞一下 1.准备需要的一些数据：1.1 从中国城市温室气体工作组下载相应的数据： 笔者存储在本地为2005.xlsx等，如下所示：  1.2 中国地级市行政区划矢量数据（空间匹配）、中国省级行政区划矢量和九段线矢量（绘图）这个学地理相关的应该都有吧（没有的私信笔者俺可以发给你） 注：为了方便空间匹配，我已提前整理了一份同碳排放数据城市名对应的标准市级行政区划名称表格 2.空间匹配碳排放数据和地级市行政区划矢量数据做空间匹配以获取空间信息 导入库 library(sf) library(readxl) library(writexl) library(tidyverse)读取地级市行政区划矢量数据和整理好的标准城市名表格数据（此表格和下载到的碳排放数据的城市一一对应） city= read_sf('./data/city/') |> st_transform(4326) address = read_xlsx('./data/address.xlsx')cityghg提供的城市名不知道是啥时候的，还有河南、湖北和海南直辖县这三个城市，这些城市在我拥有的地级市行政区划矢量里没有，我们需要将相应的区域合成： address %>% filter(str_detect(city,'直辖')) 河南直辖的是济源市我已提前在城市名表格（即addresss.xlsx中处理成济源市） # 处理湖北 city %>% filter(省=='湖北省') %>% filter(类型=='省直辖县') 碳排放数据中的湖北直辖县应该是仙桃、潜江和天门 整理一下： # 处理湖北 city %>% filter(省=='湖北省') %>% filter(类型=='省直辖县') %>% filter(市!='神农架林区') %>% st_union() %>% st_sf(省='湖北', 市='湖北直辖', 类型='省直辖县', geometry=.) -> hubei.zx city %>% filter(省=='湖北省') %>% filter(类型!='省直辖县') %>% filter(市!='神农架林区') %>% select(省,市,类型) %>% rbind(hubei.zx) -> hubei 对于海南，三沙市包括了好多岛屿，笔者是直接把海南的省直辖县合并到了一起： # 处理海南 city %>% filter(省=='海南省') %>% filter(类型=='省直辖县') %>% st_union() %>% st_sf(省='海南', 市='海南直辖', 类型='省直辖县', geometry=.) -> hainan.zx city %>% filter(省=='海南省') %>% filter(类型!='省直辖县') %>% filter(市!='三沙市') %>% select(省,市,类型) %>% rbind(hainan.zx) -> hainan 然后将地级市的矢量面重新组合： # 矢量面重新组合： city %>% filter(!省 %in% c('湖北省','海南省')) %>% select(省,市,类型) %>% rbind(hubei,hainan) -> carbon.polygon  连接地址和矢量面 # 连接地址和矢量面 address %>% left_join(carbon.polygon, by=c('city'='市')) %>% select(city,province,geometry) -> address设置为sf类（投影为阿尔伯斯投影） albers = "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105" st_set_geometry(address,address$geometry) %>% st_transform(st_crs(albers)) -> polygon计算相应面的质心和面积： polygon$area = units::set_units(st_area(polygon), km^2) st_geometry(polygon) %>% st_centroid() %>% st_transform(4326) %>% st_coordinates() -> polygon[c('long','lat')]查看一下匹配到的矢量面： polygon plot(polygon[1])  3.绘图3.1省级行政区划和九段线底图数据读取（并转换投影）：province = read_sf('./data/province.geojson') |>
st_transform(st_crs(albers))
nineline = read_sf('./data/nineline.geojson') |>
st_transform(st_crs(albers))
province = st_union(province,nineline)3.2 原始数据为2005、2010、2015、2020四年，直接定义个绘图函数方便些：library(ggspatial)# 绘制比例尺、指北针
mapping = function(t){
filename = paste0('./data/carbon/',t,'.xlsx')
polygon$carbon=read_xlsx(filename,
sheet = '二氧化碳排放(万吨)')$总排放Total
fig1 = ggplot() +
geom_sf(data = polygon,
mapping = aes(fill=carbon,color=carbon))+
geom_sf(size = .2, fill = "transparent", color = "#060d1b", data=province)+
scale_fill_gradientn(colours = colors)+
scale_color_gradientn(colours = colors)+
labs(fill = "Carbon", color = "Carbon")
fig2 = fig1+
coord_sf(crs = st_crs('epsg:4326')) + ## 将投影坐标转换为大地坐标
scale_x_continuous(expand = c(0, 0), limits = c(107, 122), breaks = seq(70, 140, 10)) +
scale_y_continuous(expand = c(0, 0), limits = c(2, 24), breaks = seq(10, 60, 10)) +
guides(fill = "none", color = "none") +
theme_bw() +
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank()
)
fig1 +
coord_sf(crs = st_crs(albers), default_crs = st_crs('epsg:4326')) +
scale_x_continuous(expand = c(0, 0),
limits=c(72,142),
breaks=seq(70, 140, 10)) +
scale_y_continuous(expand = c(0, 0),
limits = c(17,55.5),
breaks = seq(10, 60, 10)) +
annotation_scale(location = "bl") +
annotation_north_arrow(location = "tl", style = north_arrow_nautical(
fill = c("grey40", "white"), line_col = "grey20")) +
annotation_custom(ggplotGrob(fig2),
xmin= 122,xmax = 138,
ymin=15,ymax = 29)+
theme_bw() +
theme(
panel.grid=element_blank(),
legend.position = c(0.99,0.85),
legend.justification = c(1,1),
axis.text= element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank()
)-> fig
return(fig)
}定义使用的总碳排放量可视化的色带： colors=c( '#2166ac', '#99d594', '#e6f598', '#fee08b', '#fc8d59', '#d53e4f' )3.3 绘制四年的数据并用patchwork组合成一张图：f1 = mapping(2005) f2 = mapping(2010) f3 = mapping(2015) f4 = mapping(2020)组合成图： library(patchwork) f1 + labs(title="(a)2005")+ f2 + labs(title="(b)2010")+ f3 + labs(title="(c)2015")+ f4 + labs(title="(d)2020")+ plot_layout(nrow=2, widths = c(8,8), heights = c(6,6))-> fig 保存图片： ggsave("./figure/carbon2020.png", fig, height=12,width=16, dpi=600)总的代码如下： library(sf) library(readxl) library(writexl) library(tidyverse) city= read_sf('./data/city/') |> st_transform(4326) address = read_xlsx('./data/address.xlsx') address %>% filter(str_detect(city,'直辖')) # 处理湖北 city %>% filter(省=='湖北省') %>% filter(类型=='省直辖县') %>% filter(市!='神农架林区') %>% st_union() %>% st_sf(省='湖北', 市='湖北直辖', 类型='省直辖县', geometry=.) -> hubei.zx city %>% filter(省=='湖北省') %>% filter(类型!='省直辖县') %>% filter(市!='神农架林区') %>% select(省,市,类型) %>% rbind(hubei.zx) -> hubei # 处理海南 city %>% filter(省=='海南省') %>% filter(类型=='省直辖县') %>% st_union() %>% st_sf(省='海南', 市='海南直辖', 类型='省直辖县', geometry=.) -> hainan.zx city %>% filter(省=='海南省') %>% filter(类型!='省直辖县') %>% filter(市!='三沙市') %>% select(省,市,类型) %>% rbind(hainan.zx) -> hainan # 矢量面重新组合： city %>% filter(!省 %in% c('湖北省','海南省')) %>% select(省,市,类型) %>% rbind(hubei,hainan) -> carbon.polygon # 连接地址和矢量面 address %>% left_join(carbon.polygon, by=c('city'='市')) %>% select(city,province,geometry) -> address albers = "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105" st_set_geometry(address,address$geometry) %>% st_transform(st_crs(albers)) -> polygon # 面积:平方千米 polygon$area = units::set_units(st_area(polygon), km^2) st_geometry(polygon) %>% st_centroid() %>% st_transform(4326) %>% st_coordinates() -> polygon[c('long','lat')] library(ggspatial) province = read_sf('./data/province.geojson') |> st_transform(st_crs(albers)) nineline = read_sf('./data/nineline.geojson') |> st_transform(st_crs(albers)) province = st_union(province,nineline) colors=c( '#2166ac', '#99d594', '#e6f598', '#fee08b', '#fc8d59', '#d53e4f' ) mapping = function(t){ filename = paste0('./data/carbon/',t,'.xlsx') polygon$carbon=read_xlsx(filename, sheet = '二氧化碳排放(万吨)')$总排放Total fig1 = ggplot() + geom_sf(data = polygon, mapping = aes(fill=carbon,color=carbon))+ geom_sf(size = .2, fill = "transparent", color = "#060d1b", data=province)+ scale_fill_gradientn(colours = colors)+ scale_color_gradientn(colours = colors)+ labs(fill = "Carbon", color = "Carbon") fig2 = fig1+ coord_sf(crs = st_crs('epsg:4326')) + ## 将投影坐标转换为大地坐标 scale_x_continuous(expand = c(0, 0), limits = c(107, 122), breaks = seq(70, 140, 10)) + scale_y_continuous(expand = c(0, 0), limits = c(2, 24), breaks = seq(10, 60, 10)) + guides(fill = "none", color = "none") + theme_bw() + theme( axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank() ) fig1 + coord_sf(crs = st_crs(albers), default_crs = st_crs('epsg:4326')) + scale_x_continuous(expand = c(0, 0), limits=c(72,142), breaks=seq(70, 140, 10)) + scale_y_continuous(expand = c(0, 0), limits = c(17,55.5), breaks = seq(10, 60, 10)) + annotation_scale(location = "bl") + annotation_north_arrow(location = "tl", style = north_arrow_nautical( fill = c("grey40", "white"), line_col = "grey20")) + annotation_custom(ggplotGrob(fig2), xmin= 122,xmax = 138, ymin=15,ymax = 29)+ theme_bw() + theme( panel.grid=element_blank(), legend.position = c(0.99,0.85), legend.justification = c(1,1), axis.text= element_blank(), axis.ticks = element_blank(), axis.title = element_blank() )-> fig return(fig) } f1 = mapping(2005) f2 = mapping(2010) f3 = mapping(2015) f4 = mapping(2020) library(patchwork) f1 + labs(title="(a)2005")+ f2 + labs(title="(b)2010")+ f3 + labs(title="(c)2015")+ f4 + labs(title="(d)2020")+ plot_layout(nrow=2, widths = c(8,8), heights = c(6,6))-> fig ggsave("./figure/carbon2020.png", fig, height=12,width=16, dpi=600)完结收工！ |
【本文地址】