如何利用shp文件和遥感 土地利用地图进行落图

华东师范大学地图学与地理信息系统(GIS)考研试题____百度文库
两大类热门资源免费畅读
续费一年阅读会员,立省24元!
华东师范大学地图学与地理信息系统(GIS)考研试题___
上传于||暂无简介
阅读已结束,如果下载本文需要使用2下载券
想免费下载本文?
下载文档到电脑,查找使用更方便
还剩4页未读,继续阅读
你可能喜欢基于高分辨率遥感影像的农村聚落信息的提取_图文_百度文库
两大类热门资源免费畅读
续费一年阅读会员,立省24元!
基于高分辨率遥感影像的农村聚落信息的提取
上传于||文档简介
&&基​于​高​分​辨​率​遥​感​影​像​的​农​村​聚​落​信​息​的​提​取
阅读已结束,如果下载本文需要使用2下载券
想免费下载本文?
下载文档到电脑,查找使用更方便
还剩3页未读,继续阅读
你可能喜欢当前位置: &
2,315 次阅读 -
【注】新版本的maptools包对很多函数进行了修改,对于修改的内容,文章中用红色的文字进行了说明。
鉴于最近有不少人在讨论用R软件绘制地图的问题,我也就跟着凑了凑热闹,对相应的方法学习了一番。下面的这篇文章是一个初步的介绍,还有很多内容仍在学习和探索中,如果大家有什么意见或建议,我将根据自己学习的情况对文章进行进一步的补充。
在R中绘制地图其实是十分方便的,最直接的办法大概就是安装maps和mapdata这两个包,然后输入下面的命令:
library(maps)
library(mapdata)
map("china")
library(maps)
library(mapdata)
map("china")
其中map()函数还可以加上很多参数,在这里就不一一详述,具体的用法只需问号之。然而仔细看一看这张地图你会发现重庆市和四川省仍然是浑然一体,可见该地图的数据应该是有些年头了。
幸运的是,通过谢益辉的这篇博文我们已经可以大体知道该如何操作了,下面就为大家介绍一下具体的步骤。
首先,从这里下载中国地图的GIS数据,这是一个压缩包,完全解压后包含三个文件(bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx),将这三个文件解压到同一个目录下,并在R中设好相应的工作空间,然后安装maptools包,运行如下程序:
library(maptools);
x=read.shape('bou2_4p.shp');#下文中会继续用到x这个变量,
#如果你用的是其它的名称,
#请在下文的程序中也进行相应的改动。
library(maptools);
x=read.shape('bou2_4p.shp');#下文中会继续用到x这个变量,&&
&&&&&&&&&&&&&&&&&&&&&&&& #如果你用的是其它的名称,&&
&&&&&&&&&&&&&&&&&&&&&&&&&&#请在下文的程序中也进行相应的改动。
【修改】新版本的maptools包不再提供read.shape()函数,请用readShapePoly()代替。
这时一张完整的中国地图就已经画好了。但是在实际使用的过程中,我们往往会根据自己的需要对地图中的某些省份着以特定的颜色,这时就可以通过调节plot命令中的fg参数来予以实现。然而为了清楚地说明这部分的内容,我需要插播一段R绘制地图的原理。
在绘制地图时,每一个省市自治区或者岛屿都是用一个多边形来表示的。之前的GIS数据,其实就是提供了每一个行政区其多边形逐点的坐标,然后R软件通过顺次连接这些坐标,就绘制出了一个多边形区域。在上面的数据中,一共包含了925个多边形的信息,之所以有这么多是因为一些省份有很多小的附属岛屿。在这925个多边形中,每一个都对应一个唯一的ID,编号分别从1到925。
回到刚才的话题,plot命令中的fg参数在本例中应该是一个长度为925的向量,其第i个分量的取值就代表了地图中第i个多边形的颜色。一个简单的尝试是运行下面这个命令看看效果:
plot(x,fg=gray(924:0/924));
plot(x,fg=gray(924:0/924));
【修改】新版本的maptools包的绘图参数也有所改变,请将fg换成col。
于是自然就产生了一个问题:如何获取某一个特定地区的ID,进而设置我们想要的颜色?事实上,在变量x中,就已经存储了我们想要的信息。在R中输入“x[[2]]”或“x$att.data”,会得到一个925行7列的数据框,这其实是bou2_4p.dbf这个文件中存储的信息,之前的read.shape()函数虽然读取的是bou2_4p.shp文件,但在默认情况下会把dbf文件的信息也放到变量之中。对于这个数据框,其行名就是每一个区域的ID编号,第一列和第二列分别是面积和周长,最后一列是该区域所属的行政区名,其它的列应该也是一些编号性质的变量。于是,通过查找相应的行政区对应的行名,就可以对fg参数进行赋值了。下面是我编的一个函数,用来生成所需的fg向量:
getColor=function(mapdata,provname,provcol,othercol)
f=function(x,y) ifelse(x %in% y,which(y==x),0);
colIndex=sapply(mapdata$att.data$NAME,f,provname);
fg=c(othercol,provcol)[colIndex+1];
return(fg);
getColor=function(mapdata,provname,provcol,othercol)
f=function(x,y) ifelse(x %in% y,which(y==x),0);
colIndex=sapply(mapdata$att.data$NAME,f,provname);
fg=c(othercol,provcol)[colIndex+1];
return(fg);
【修改】地图数据的组织形式有所变化,上面函数中的mapdata$att.data$NAME需要替换为mapdata@data$NAME。
其中mapdata是存放地图数据的变量,在上面的例子中就是x,provname是需要改变颜色的地区的名称,provcol是对应于provname的代表颜色的向量(名称和数字均可),othercol是其它地区的颜色。举例如下:
provname=c("北京市","天津市","上海市","重庆市");
provcol=c("red","green","yellow","purple");
plot(x,fg=getColor(x,provname,provcol,"white"));
provname=c("北京市","天津市","上海市","重庆市");
provcol=c("red","green","yellow","purple");
plot(x,fg=getColor(x,provname,provcol,"white"));
注意provname一定要写地区的全称,写法可以参照下面这条命令生成的向量:
as.character(na.omit(unique(x$att.data$NAME)));
as.character(na.omit(unique(x$att.data$NAME)));
由此生成的向量有33个元素,少了澳门特别行政区,这是这个数据中的一块瑕疵。在x$att.data的第899行有一个NA,不知道它代表的是否就是澳门。
利用类似的方法就可以根据自己的需要对不同的区域进行着色,下面再举一例。从国家统计局获取2007年我国各地区的人口数据,然后根据人口的多少对各省份进行着色。程序如下:
provname=c("北京市","天津市","河北省","山西省","内蒙古自治区",
"辽宁省","吉林省","黑龙江省","上海市","江苏省",
"浙江省","安徽省","福建省","江西省","山东省",
"河南省","湖北省","湖南省","广东省",
"广西壮族自治区","海南省","重庆市","四川省","贵州省",
"云南省","西藏自治区","陕西省","甘肃省","青海省",
"宁夏回族自治区","新疆维吾尔自治区","台湾省",
"香港特别行政区");
pop=c(43,98,58,7625,
6,14,284,,
552,610,3);
provcol=rgb(red=1-pop/max(pop)/2,green=1-pop/max(pop)/2,blue=0);
plot(x,fg=getColor(x,provname,provcol,"white"),xlab="",ylab="");
1234567891011121314
provname=c("北京市","天津市","河北省","山西省","内蒙古自治区",
"辽宁省","吉林省","黑龙江省","上海市","江苏省",
"浙江省","安徽省","福建省","江西省","山东省",
"河南省","湖北省","湖南省","广东省",
"广西壮族自治区","海南省","重庆市","四川省","贵州省",
"云南省","西藏自治区","陕西省","甘肃省","青海省",
"宁夏回族自治区","新疆维吾尔自治区","台湾省",
"香港特别行政区");
pop=c(1633,1115,6943,3393,2405,4298,2730,3824,1858,7625,
5060,6118,3581,4368,9367,9360,5699,6355,9449,
4768,845,2816,8127,3762,4514,284,3748,2617,
552,610,2095,2296,693);
provcol=rgb(red=1-pop/max(pop)/2,green=1-pop/max(pop)/2,blue=0);
plot(x,fg=getColor(x,provname,provcol,"white"),xlab="",ylab="");
其中颜色越深的地方代表人口数越多,反之为人口数越少。
此外,在绘制地图的过程中,还有一个比较有用的参数是recs,它是一个由多边形ID组成的向量,表示在地图中只画出这些ID所代表的区域。利用这个参数,就可以画出某一部分的地图,例如下面的例子是我国中部六省的地图:
getID=function(mapdata,provname)
index=mapdata$att.data$NAME %in%
ids=rownames(mapdata$att.data[index,]);
return(as.numeric(ids));
midchina=c("河南省","山西省","湖北省","安徽省","湖南省","江西省");
plot(x,recs=getID(x,midchina),fg="green",ol="white",xlab="",ylab="");
getID=function(mapdata,provname)
index=mapdata$att.data$NAME %in% provname;
ids=rownames(mapdata$att.data[index,]);
return(as.numeric(ids));
midchina=c("河南省","山西省","湖北省","安徽省","湖南省","江西省");
plot(x,recs=getID(x,midchina),fg="green",ol="white",xlab="",ylab="");
上面的getID()是我编写的一个功能与getColor()类似的函数,用来返回指定省份的ID。
【修改】新版本的maptools包的绘图函数已经取消了recs这个参数,现在要实现这个功能,可以在颜色上把不需要的省份变成白色,其中填充色用col参数,边界颜色用border参数。例如上面的例子可以用下面的函数来实现:
plot(x, col = getColor(x, midchina, rep("green", 6),
"white"), border = "white", xlab = "", ylab = "")
plot(x, col = getColor(x, midchina, rep("green", 6),
&&&&"white"), border = "white", xlab = "", ylab = "")
最后要说的是,在画出的图上仍然可以用points()函数和text()函数加上点和文字,而maptools包中还提供了一个pointLabel()函数,用来解决文本标签的重叠问题。这部分内容请参阅博文:用R画中国地图并标注城市位置,以及避免文本标签重叠:maptools中的pointLabel()。
从以上的内容来看,本文所述的都是一些最基本的绘图方法,还没有对地理信息数据进行更进一步的分析。如果有机会的话,这一主题的下一篇文章将为大家介绍地图数据的组成结构,并说明如何将不同格式的地理数据整合起来,例如如何在上面的地图上绘制出我国的铁路、水系分布等内容。
原文链接:/article/0I133X2014.html
注:转载文章均来自于公开网络,仅供学习使用,不会用于任何商业用途,如果侵犯到原作者的权益,请您与我们联系删除或者授权事宜,联系邮箱:contact@dataunion.org。转载数盟网站文章请注明原文章作者,否则产生的任何版权纠纷与数盟无关。
相关文章!
不用想啦,马上 发表自已的想法.
做最棒的数据科学社区
扫描二维码,加微信公众号
联系我们:没想到上个星期在写一篇论文(非毕业论文,)的时候,居然需要画地图…这个,地图怎么画?虽然常年看论文人家总是轻描淡写的说一句“GIS数据来自……”,我还是对这个东西没啥概念,总觉得貌似挺麻烦似的……
一直知道R能画,taiyun展示过,,,但是毕竟自己没有亲手做过。这次逼上梁山了,不得不搞定。于是乎,照着葫芦画瓢,比着yixuan和yihui的教程,一步步的研究,研究。
因为这次我画地图主要是画中国地图,所以GIS数据自然是关注中国的。从yihui那里找到了“”(中文版和我一样不给力的同学请用英文版网址:)。然后,嗯,GIS数据就有了。扔到R里面,调用maptools包,两行代码,中国的雄鸡图就出来啦!顿时那个爽啊。
可是然后呢?呃,我要标注的是我所感兴趣的几个城市的位置,这个yixuan在“可能的拓展”里面说了一下,没有说具体的(傻瓜式)操作,yihui的地图用的是文本导入经纬度数据,我可没有现成的经纬度数据啊!我又懒得研究这个GIS对象在R中到底是怎么玩的,就对着maptools包里面的各种函数扫过去。果然,有一个函数貌似是可以读取经纬度数据的。联系起国家基础地理信息系统里面是有提供各个城市的经纬度数据的,嗯,试了试这个函数,果然可以!剩下的工作就轻松了,几行程序告诉R那几个城市是我感兴趣的,然后一个命令就上去了,再一个命令就搞定名称的一并标注了。R果然聪明啊!
成品图在此:
好了,最后简单的说一下我用到的代码,方便大家复制粘帖。
library(maptools); #调用maptools包
#read data
cities=readShapePoints('chinamap/res2_4m.shp') #国家地理信息系统下载的市级经纬度数据
china=readShapePoly('chinamap/bou2_4p.shp') #国家地理信息系统下载的省级多边形数据
plot(china); #画中国地图
points(cities[cities$NAME %in% piaohao_pre1850$city,], pch = 19, col = rgb(0, 0, 0, 0.5)) #标注感兴趣的城市黑点
text(cities[cities$NAME %in% piaohao_pre1850$city,],labels=cities[cities$NAME %in% piaohao_pre1850$city,]$NAME, cex = 0.9, col = rgb(0,
0, 0, 0.7))
#标注城市名称
顺带赞一个,国家地理信息系统里面城市名分别使用中文和拼音存储的,调用的时候任意一个都可以。我的代码里面piaohao_pre1850$city存的是我感兴趣的城市的名称,只要对应一下就可以了,很方便。
(另注:汇通天下,嘿嘿~)
这票号数据哪里找的啊?
一点点声明
怕下次找不到门?直接google“落园”呗。
落园是我的非学术博客,只是为了娱乐大众。如果您对学术感兴趣,请移步我的或查看我的。
长篇连载系列

我要回帖

更多关于 土地利用动态遥感调查 的文章

 

随机推荐