||
平常的中国地图,南海部分是一个独立的小图,放在整张地图的一角的。不知道专业的地图软件是不是很容易实现这样的效果,反正我是费了整整一个工作日的时间才在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)) #本句与主题无关。让做图参数退回默认值。
#-----------程序结束
只要你有中国边界的数据,可以把我的代码直接拷贝执行。也可以把注释都删除再执行。
我的编辑器和输入法真是不给力(老吞键),把这段文字打完,累出一身汗,愣是把脑力劳动变成体力劳动了。虽然言猶未尽,也只好就此罢手。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-24 01:44
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社