0. 今日绘图目标
将map与柱状图绘制在一张图中。
1. 所需R-packages
require(ggplot2)
require(raster)
require(ncdf4)
require(reshape2)
#软件包安装方式
install.packages('ggplot2')
install.packages('raster')
install.packages('ncdf4')
install.packages('reshape2')
2. 示例数据
2.1 ret_box #包括long, lat, cor(用于控制空间点的大小)
long lat cor
1 319.9667 -21.88330 -0.2909924846
2 233.6058 49.80639 -0.6612940908
4 322.1500 -33.38330 -0.6554984115
5 320.2167 -34.81670 -0.6140925293
6 322.2167 -36.71670 -0.6682152245
7 289.8667 46.26670 -0.2671691399
8 67.7900 34.12000 -0.4361970541
9 74.4850 32.02500 -0.6861450533
10 74.5330 38.35000 -0.8576893172
11 60.0500 42.05000 -0.7424349943
12 74.4500 38.26700 -0.8279314389
13 66.7830 39.70700 -0.7111693976
14 253.0167 32.65000 -0.6827918702
15 296.3333 48.80000 0.4564051335
16 253.4333 25.76670 -0.7326677653
17 271.9333 48.36670 -0.2170552270
18 251.5667 32.96670 -0.6881968791
19 113.1503 -17.31530 0.7252337143
20 113.6964 -20.96580 0.8489313090
21 75.1170 31.43300 -0.7266186825
22 63.8330 39.78300 -0.6851615075
23 66.3940 38.33600 -0.7178798947
24 61.3000 48.05000 0.4712528403
25 70.0900 31.43000 -0.7278034435
26 66.9370 38.45000 -0.7069809031
27 75.6700 29.81000 -0.5334788720
28 67.6990 45.22400 -0.5065231875
29 67.5750 45.17200 -0.7963865354
30 67.6910 45.22800 -0.5065231875
31 62.4830 36.53300 0.5342509122
32 63.2670 40.63300 -0.8255034957
33 63.0990 40.28100 -0.7726425553
34 62.9590 40.31500 -0.7262693955
35 63.2810 40.59900 -0.8237684070
36 61.7810 38.01100 0.2010569215
38 61.2030 40.26700 0.0100699356
40 271.7167 50.41670 -0.8218614332
41 273.3333 50.33330 -0.6497722055
42 251.3333 33.33330 -0.4800369132
43 252.7500 32.60000 -0.6543348992
45 254.3569 47.53861 -0.2109533665
46 70.5090 40.60700 -0.4188613454
47 63.5670 36.45000 0.5692960424
48 71.8300 44.58000 0.1195021740
49 71.8300 44.58000 0.1195021740
50 63.7500 35.95000 0.3014936360
2.2 全球地图数据－格式 .shp
全球地图下载地址
注：该数据直接导入R是不显示坐标的
class : SpatialPolygonsDataFrame
features : 7
extent : -180, 180, -55.72057, 83.62303 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 1
names : Region2
min values : Africa
max values : South America
3. 数据预处理
3.1 ret_box 预处理
ret_m = melt(ret_box,c('long','lat'))
ret_m$long = ret_m$long - 180
index_neg = which(ret_m$value <=0)
index_pos = which(ret_m$value >0)
ret_pos = ret_m[index_pos,]
ret_neg = ret_m[index_neg,]
ret_pos$cuts = cut(ret_pos$value, breaks = c(0,0.25,0.50,0.75,1.00))
ret_neg$cuts = cut(ret_neg$value, breaks = c(-1.00,-0.75,-0.5,-0.25,0))
shape_neg = rep(25,length(index_neg))
shape_pos = rep(24,length(index_pos))
ret_pos$shape = shape_pos
ret_neg$shape = shape_neg
ret_box = rbind(ret_neg,ret_pos)
ret_box$type = c(rep(2,nrow(ret_neg)),rep(3,nrow(ret_pos)))
ret_box2 = ret_box
ret_box2$type = 1 #控制柱状图在第一个facet 中
neg_index = which(ret_box2$value <0)
pos_index = which(ret_box2$value >0)
ret_box2$fill = 'blue'
ret_box2$fill[neg_index] = 'red'
3.1 Global map数据预处理
world = shapefile('...地图路径地址...')
prov_map = fortify(world,region = "Region2") #Region2 为地图中类的名称
prov_map2 = rbind(prov_map,prov_map)
prov_map2$type = rep(c(2,3),each = nrow(prov_map)) # 2，3 为画在第几个分面中
4. 正式绘图
color_map = c('#E8002D','#FE2725','#FE5B34','#FEB666',
'#6CBDE3','#177ABD','#2A4FA0','#013F87')
a = ggplot()+
geom_histogram(data = ret_box2, aes(x = value,fill = fill),binwidth = 0.25,center = 0.125)+
geom_polygon(data=prov_map2,aes(x=long, y= lat,group = group),
colour='black',fill = 'transparent',size = 0.5 )+
coord_fixed(1.3)+
geom_point(data = ret_box,aes(x = long, y = lat,color = cuts,size = cuts))+
scale_color_manual(values = color_map)+
facet_wrap(~type, nrow = 3,scales = 'free')
png('...图件输出路径...',
width = 20,
height = 20,
units = 'cm',
res = 1000)
print(a)
dev.off()
5. 出图示例
最终图件示例
6. 技术支持-微信
大家有看不懂的地方欢迎在文下评论，或者可以私戳我哈~
微信联系方式
网友评论