本帖最后由 R语言微信号 于 2015-3-10 18:12 编辑
2015-03-04xxiao513
朋友圈分享 并回复：0304 即可获取本节程序和R数据文件
上期内容回顾上一期中主要介绍了地理信息系统(GIS)及其所处理的空间数据类型。
温馨提示：如果您对此感到有些陌生，那就赶快翻开公众号历史消息温习一下，点击文章末尾“阅读原文”可见~温故而知新哦~ 本期内容导读本期主要介绍空间数据如何在R语言的中sp对象和GIS通常所使用的外部格式数据对象间进行转换。首先介绍坐标系的转换，在此基础上介绍如何快捷地使用坐标参考系进行空间数据的导入和导出。 ■ 1.坐标参考系（CRS）和投影 坐标参考系是记录空间位置信息的基础，空间数据常用的坐标系分为两大类：地理坐标系（Geographiccoordinate system，为球面坐标，参考平面为椭球面，坐标单位为经纬度）和投影坐标系（ProjectedCoordinate Systems，为平面坐标，参考平面为水平面，坐标单位为米、千米等）。通常我们看到的地图是将地理坐标转换到投影坐标的结果，这个转换过程称为投影。 在R中，多种不同坐标系下的数据如何通过坐标系转换从而得到统一的表达和可视化非常重要。这一步用到的核心包是sp和rgdal，分别用于定义、读取并转换空间坐标系。在R中CRS由PROJ.4 库中的proj4string即一个字符串定义，以中国省级地图为例我们看一个地图坐标转换和投影的例子： #加载包 library(rgdal) library(sp) #读入中国省级地图,为一个SpatialPolygonsDataFrame对象 China <- readOGR("D:\\China_province_shp", "Provinces_region")#获取CRS信息 proj4string(China)[1] "+proj=longlat +ellps=WGS84 +no_defs"
CRS("+init=epsg:4030")CRS arguments:
+init=epsg:4030 +proj=longlat +ellps=WGS84 +no_defs
可以看到参数和所获取的一模一样，这个坐标系是一个以经纬度为单位的球面地理坐标系，如果以经纬度为横纵坐标绘图，绘图结果如下： plot(China, axes=TRUE, border="gray") 接下来，将球面坐标投影成平面坐标，并添加经纬度格网投影的结果： #定义一个新的平面投影坐标系，叫北京54 newCRS <- CRS("+init=epsg:2431")#将地图坐标进行投影变换得到新地图 China_Proj <- spTransform(China, newCRS)#查看投影信息 proj4string(China_Proj)#可视化 plot(China_Proj, axes=TRUE, border="grey")#获取地理坐标系下的经纬度格网 grd_LL <- gridlines(China, ndiscr=100)#将格网投影为二维平面坐标 grd_BNG <- spTransform(grd_LL, newCRS)#添加投影后的格网 plot(grd_BNG, add=TRUE, lty=2)#获取球坐标系下网格线的经纬度标记点 grdtxt_LL <- gridat(China)#将网格线经纬度标记转换为投影坐标系下 grdtxt_BNG <- spTransform(grdtxt_LL, newCRS)#添加投影后的标记 text(coordinates(grdtxt_BNG), labels=parse(text=as.character(grdtxt_BNG$labels)))绘图结果如下所示，可以看到球坐标到平面坐标存在着变形，所以在可视化时需要选择好坐标系并将数据统一到统一坐标系下。■ 2.空间数据的导入和导出 第一期我们提到空间数据分成两大类，一类是矢量数据包括点线面等，第二类是栅格数据，如卫星地图影像。 2.1 矢量数据导入导出 常使用的矢量空间数据格式是ESRI规定的shapefile。这种格式用至少三个文件表示矢量数据，其中最主要的就是存储几何拓扑信息的.shp文件和存储属性数据的.dbf文件。读入该数据格式的方法有多种，在1部分就介绍了通过rgdal包的readOGR()函数读入外部shp格式的中国地图（世界各国的行政区划文件可以从这里http://gadm.org/找到）。 library(rgdal) #读入中国省级地图,为一个SpatialPolygonsDataFrame对象 China <- readOGR("D:\\China_province_shp","Provinces_region") 除rgdal包外，还有maptools包的readShapePoints、readShapelines和readShapePoly函数，用于读取shp数据到R的sp对象。另外还有shapefile包、PBSmapping包等也可以以类似的方式读入shp数据转换为R的空间对象。 导出和一般的R文件导出一样，有read就有对应的write。例如我们想将1部分中投影后的中国地图保存为shp文件： #输出投影后的中国地图为shp文件 writeOGR(China_Proj, "D:\\China_province_shp", "Provinces_region_prj", driver="ESRI Shapefile")也可以将地理坐标系下的中国省级地图输出为谷歌地球能加载的kml格式，这样就能以更有趣的互动方式查看数据了。 #输出地理坐标系下的中国地图为kml文件 writeOGR(China, "D:\\China_province_shp\\ChinaProvinces.kml", "Provinces_region", driver="KML")在安装完谷歌地球后打开，加载输出的kml文件：
2.2 栅格数据导入导出 栅格数据是图像数据，其中每个栅格内有其属性值，最常见的栅格数据为高程数据，记录每个位置的海拔高度值。处理栅格数据的包主要有raster和rgdal包，raster的读取和输出部分依赖于rgdal包。以Applied Spatial Data Analysis with R书中的奥克兰海拔数据为例导入格式为tiff的栅格数据并可视化： library(raster)#使用raster函数读入tif格式的栅格数据 r <- raster("D:\\China_province_shp\\70042108.tif")#将栅格数据中海拔小于等于0的部分值设置为NA，绘图不显示 out <- raster(r)bs <- blockSize(out)out <- writeStart(out, filename = tempfile(), overwrite = TRUE)for (i in 1:bs$n) { v <- getValues(r, row = bs$row, nrows = bs$nrows) v[v <= 0] <- NA writeValues(out, v, bs$row) }out <- writeStop(out)par(mar=c(2,0.5,0.5,1))plot(out, col = terrain.colors(100))将去除海平面以下高度数据的海拔数据输出为tif格式：rf <- writeRaster(out, filename="D:\\China_province_shp\\test.tif", format="GTiff", overwrite=TRUE)
下一期预告 将为大家介绍空间数据可视化部分。先从传统的R语言绘图系统开始，介绍更加灵活的图形可视化技巧。
关注我们—官方网站— —官方QQ群— R语言中文论坛-2（1000人群）：427060123 R语言中文论坛（2000人群，已满）：74076289 Biostatistician（500）：186701945 —官方微博— 新浪微博：@R语言中文网官网 —官方微信— 微信名：R语言中文网 微信号：rchinanet
点击“阅读原文”温习上期内容！喜欢就点个赞鼓励一下~
关注该公众号
| 
发表于 2015-3-10 18:10:37