牛海山博客分享 http://blog.sciencenet.cn/u/leymus

博文

用R画中国地图时如何把南海做成小图

已有 9328 次阅读 2014-12-1 14:08 |个人分类:R使用|系统分类:科研笔记

平常的中国地图,南海部分是一个独立的小图,放在整张地图的一角的。不知道专业的地图软件是不是很容易实现这样的效果,反正我是费了整整一个工作日的时间才在R里做到这一点。赶紧上来宣传一下,一个是表表功,再一个是记录下来、留作后用。


其实,我以前想都没想过为什么要把中国地图折成这样的两截。也许审美观随着年龄增长吧,我今天突然产生这样的需求了。因为我想表现中国陆地上的某一种趋势,如果把南海部分疆土放进去吧,整个地图中陆地比例就太小了、不好看,如果把南海部分切下去吧,又于心不忍。于是,就想起来用传统的“两截图”了。


在R里如何实现呢?上网查了很久,也没有找到答案。有人用GrADS很容易就实现这种效果了,寥寥几行命令,看着都爽。心中暗想:将来一定要学学这个简介明快的语言。但是远水不解近渴,暂时还不能用它。也有传授用R画中国地图的,要么直接把南海略去不表了,要么画个竖版地图,——这样的方案不能接受(叹号)。第三条道路,就是采用 http://www.gadm.org/  提供的中国行政边界图,这个图里面就没有南海部分。这个数据库也是权威发布,很多作者都采用。本来我用一用也是无妨的,但在好奇心的召唤下,我还是看了看这个版本的中国地图除了把南海切下去了之外,与国家测绘局的“标准”中国地图还有什么区别。这一比之下,才发现差别不小。首先,是麦克马洪线南边(也就是中印争端东段)——比台湾省面积还大的——那部分领土都不算中国的了。它按实际控制区划的线。其次,这个图与中国的”标准图“还有很多边界不一致的地方,虽然很细微,因为中国很大,在地图上看不出来,但是具体数字上能显示出来。这些差异小差异也许是测量误差造成的。第三条道路,我也决定不走了。还是老老实实用R做两截图吧。
 
最初的想法是,也许R中做地图的包 (package) 中有相关命令,于是GOOGLE了一通可能的包和可能的命令。整整一个上午,我都在看文献,以前没学过地理信息系统,也没学过地图学,所以研读了大量乍看上去相关、扎进去细看却是无关的vignettes,还是没有头绪。虽然在此期间也学了不少有用的东西,眼前的问题仍然没有解决。R的包实在太多了、看不过来,真是包到用时方恨多。无奈之下,我想,既然聪明办法一时找不到,干脆用笨办法吧,两步走:先把南海部分的地图做出来,然后再把它导入、粘到大图的一个角上去。于是按照这个思路,我查了一些文档,看怎么着导入这个小图到大图上去最方便。在这个查找过程中,突然之间,眼角的余光扫到了google预览里的 par(fig=(...),new=TRUE) 这么一行代码,直觉告诉我,这一行就是通关秘诀。直觉是对的,我用了简单的几行命令就顺利解决了问题,当然解决得不是很完美、但是对于发表已经足够了,另外相对于前面提到的方案而言,也算是更进一步了。
 
下面把我的代码拷贝并解释如下,供广大领导干部和人民群众批判之用。#后面的都是注释。不包括注释,一共用了13行代码,也算是比较简洁吧? 这13行里面,还有几行并非核心代码。
 
#-----------程序开始
# 画中国地图

rm(list=ls())#清空所有对象,一切从新开始


library(sp)#是很多空间分析的基础包,需要先加载。

   #额外补充:如果你计算机上没有安装这些包,是加载不上的。

   #只有先安装这三个包(sp, raster 和 rgdal),再执行我的代码,才能在执行的过程中加载它们。
   #如果不知如何安装这三个包,问一下常用R的人,或者看一下基础文档。

   #操作系统是LINUX / UNIX 的,在安装这些包前,还要给操作系统安装一些依赖库,

   #才能在R里安装这些包。 至于需要装哪些库,编译的时候,在出错信息里面可以看到。

   #依稀记得,至少需要gdal库和proj.4。


library(raster)#投影时要用到这个包

library(rgdal)#读取矢量边界要用到此包
 
china.bou<-readOGR(dsn="./ChinaMap/bou1_4m/bou1_4p.shp",layer="bou1_4p") 
   #读入中国边界矢量数据。文件路径是我的文件路径,请按照实际情况修改。
   #这里的数据库,是国家测绘局提供的“国家基础地理信息系统数据"。
   #这是前些年下载的,我算是捡着大便宜了,因为现在这个库开始收费了。

   #拿国家边界做买卖,真是丧心病狂。

projection(china.bou)<-CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs") 
   #投影赋值,是WGS84投影

ext<-extent(china.bou)
   #算出来能框住中国的最小矩形。

   #能给出四个值(经纬度范围),分别是:x最小值,x最大值,y最小值,y最大值。

   #这四个值可用ext[1],ext[2],ext[3] 和 ext[4]提取出来,后面要用到
 
my.op<-par()    #这一句,与本题主旨无关。意思是要记住做图参数的默认值。
 
plot(china.bou,    #刚才读入的国界
   xlim=c(ext[1],ext[2]),    #定义大图中要显示的经度(x轴)范围,即x最小值和最大值
   ylim=c(ext[3]+13,ext[4])    #定义要显示的纬度范围,即y的最大、最小值。

                                                 #里面13这个参数可根据自己的审美观调整
                                                 #现在(ext[3]+13)的位置应该是海南岛稍微往下一点
)    #这一句执行完之后,中国地图的大图就出来了。当然南海部分不完整。
 
par(fig=c(.15,.35,.04,.35),    #定义小图的位置和大小。

                                               #注意其中的四个坐标用的是(0,1)区间的相对值。        
       new=TRUE,    #这一句,表示在前图上再加个图,而不是把前图抹去。

                               #我就是不知道这个参数才饶了大弯子。

                               #自认为对par()了如指掌了,其实尚未登堂入室呢。
       mar=c(0,0,0,0)    # margin的不要,统统去掉
)
 
plot(china.bou,    #再画中国地图,这次是显示在小图里面的。
   xlim=c(ext[1]+24,ext[2]-13),    #定义小图里要显示的经度范围。其中的+24,-13都可修改
   ylim=c(ext[3],ext[4]-30)   #定义小图中要显示的纬度范围。-30是可修改的参数。
)     #此句执行完,小图就画上去了。领土完整了。
 
box(lty="1373",col="red")
   #给小图加个框。”1373“代表框的类型,是我想到的最浪漫的类型,”red“是指框的颜色。

   #这俩参数可根据个人喜好修改。
 
on.exit(par(my.op))    #本句与主题无关。让做图参数退回默认值。
 
#-----------程序结束

只要你有中国边界的数据,可以把我的代码直接拷贝执行。也可以把注释都删除再执行。


我的编辑器和输入法真是不给力(老吞键),把这段文字打完,累出一身汗,愣是把脑力劳动变成体力劳动了。虽然言猶未尽,也只好就此罢手。




https://blog.sciencenet.cn/blog-39344-847675.html

上一篇:能把R下多图合并
下一篇:R做图时无法转换unicode字符的临时解决方案
收藏 IP: 124.171.86.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (3 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-11-24 01:44

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部