发帖
雷达卡
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))
下一期预告
将为大家介绍空间数据可视化部分。先从传统的R语言绘图系统开始，介绍更加灵活的图形可视化技巧。
京公网安备 11010802022788号
