Fork me 
	on GitHub

R语言读取空间数据以及ArcGIS中OLS工具回归结果可视化R语言版

前面已经介绍过R语言读取excel的方法了,当然读取数据来说,个人还是推荐csv或txt存储(针对小数据量)。大数据量的数据的话建议还是用数据库,此外也可以考虑data.table包读取,这个包也是个神包,后面学习完可能来谈谈。这个都是题外话,今天主要目的还是来介绍R语言读取空间数据的方法。主要是之前有同学问过读取的方法。我就顺带整理下,另外虾神今天刚发了一篇关于ArcGIS的OLS工具回归结果可视化内容,并贴出了Python版可视化的代码(文末贴链接),所以对应写个R语言版。
本文介绍的空间数据类型主要包括了三种:矢量数据(以最普遍的的shapefile为例),栅格(raster,这个格式就比较多了,不过大同小异),地理数据库(geodatabase也就是.gdb文件,Esri的数据库)。

1.矢量数据

矢量数据其实主要包括了三类:点,线和面,能读取的方式有很多种。下面列举几种。
先从点线面分别读取的方式来看,主要包括readShapePoints(读取点),readShapeLines(读取线要素),readShapePoly(读取面要素)。这几个函数都是maptools包里面的。所以第一步如果没安装的话请先安装。

install.packages('maptools')
library(maptools)

接着定位到我们所需读取数据的工作路径上,然后就可以开始读取对应的数据了。

fujian<-readShapePoints('fujian.shp')
nanhailine<-readShapeLines("linesnation.shp")
province<-readShapePoly("province.shp")

如果不需要什么其他操作,读取数据只需要填入文件名字作为传入的参数即可。这几个函数完整的参数大体差不多,主要包括下面几个。
fn——文件名,一般能读的是.shp文件,.shx文件和.dbf文件
proj4string = CRS(as.character(NA))——坐标系的CRS字符串,关于坐标系的问题,这里不详讲。其实就是一个坐标系对应一个ID,把对应ID读进去,按照对应坐标系读取,这个是遵循规范的。
一般前两个参数用得多。后面这些只介绍这三个函数共有的参数,其他参数就请参照帮助文档。
verbose = FALSE——默认为False,这个主要是在读取数据后是否返回读入要素的类型和数量。
repair=FALSE——这个参数的话,主要是考虑到.shx索引文件太大,默认False会跳过读取数据,TRUE的话,会进行内部修复,读取这类文件。
而maptools同样提供了另外一个函数readShapeSpatial,这个就可以读以上的三类要素。

fujian<-readShapeSpatial("fujian.shp")

当然除了maptools,还有其他包可以读取,事实上,maptools提供的函数读取只能传输较差分辨率的空间数据,所以更推荐的是用rgdal包的OGR驱动程序来读取。熟悉开源GIS的同学对GDAL会比较熟悉,事实上rgdal就是GDAL的R接口(当然没装还是要先装,方法同上),读取方式如下,参数也是传入文件名即可简单读取,不过这个参数可以读具体文件也可以读文件夹名。对应上面proj4string也有一个参数p4s,其他参数参照文档。

fujian<-readOGR("fujian.shp")

此外还有shapefiles包也可以进行读取。读取方式(可以读取shp和shx,shx读取结果为空间索引)如下:

fujian<-read.shp("fujian.shp")

矢量数据读取主要通过以上几种方式就可以实现。

2.栅格数据

栅格数据的话,格式还是多种多样的。这边主要提供几种不同格式的读取方法(.img文件,.tif文件,ASCII码文件和.asc文件)。栅格数据读取主要是基于rgdal包,读取方式如下,img和tif都可以通过readGDAL直接读取。

co2<-readGDAL("CO22008.img")
co2<-readGDAL("CO22008.tif")

这里面的参数我就不详细介绍了,主要解释几个个人认为比较重要的参数。有兴趣的同学可以去查询官方文档。
band——波段数,单纯栅格无所谓。做遥感影像数据处理时就会遇到需要几个波段的问题,如果缺省的话,是全部导入。
p4s——等同于上面的proj4string
type——像素深度:8bit,16bit等
除了rgdal之外,也可以通过raster包进行读取.img文件和.tif文件,这个更方便些。读取方式如下

co2<-raster("CO22008.img")
co2<-raster("CO22008.tif")

当然栅格数据还有较为普遍的以ASCII码文件存储的方式。这里也提供下如何读取ASCII码文件,这个方法是基于sp包的,所以需要先安装和载入sp包,这个包是R语言空间数据的基础包,指定了空间数据库的方法和对象。

co2<-read.asciigrid("co22008.txt")

当然ASCII码文件可能是以.asc文件存储的,只需把后缀名改成.asc即可读取。
栅格的读取大概就如上。

3.地理数据库(Geodatabase)数据读取

Geodatabase是Esri在ArcInfo8之后引入的一种全新的面向对象的空间数据模型。具体简介可以自己搜索。也就是说Geodatabase是Esri官方提出的一种数据库,没有ArcGIS是无法创建Geodatabase的。读取的话其实也相对麻烦些。目前看到的只有Esri官方给出的一个R包可以读取Geodatabase数据。R语言与ArcGIS的结合在未来将很有潜力。目前Esri已经在github上开源了部分工具,2015年全球用户大会上也秀出了R-ArcGIS的Sample工具。
具体开源地址:

https://r-arcgis.github.io/

本次用来读取Geodatabase的包就是R-ArcGIS中的一个关键包——arcgisbinding。这个包目前没有在cran上,建议下载之后离线安装。下载地址:

https://github.com/R-ArcGIS/r-bridge/releases/tag/v1.0.0.125。

这个包的官方文档可以从官网下载,也可以从下面的连接下载。

http://download.csdn.net/detail/esa_dsq/9823403

具体安装,同时安装完之后需要先确认ArcGIS的许可(要求应该是ArcGIS10.4以上的版本或者ArcGIS pro1.1以上),具体代码如下:

install.packages("G:/GIS/Esri/ArcGIS Plugin/arcgisbinding_1.0.0.125.zip", repos = NULL, type = "win.binary")
arc.check_product()

读取的方式稍微复杂些,用到了arc.open,arc.select,arc.data2sp三个函数,arc.open是打开gdb文件里的featureclass(支持的格式还包括layers等),arc.select是将打开的featureclass按照需要的字段和sql读成R语言中熟悉的数据框,arc.data2sp是将数据框转化成空间要素。使用方式如下。

china<-arc.open("china.gdb/province")
chinapop<-arc.select(china,fields = c('Pop_Rural','Pop_Urban','POPU'))
chinapopsp<-arc.data2sp(chinapop)

当然读取完了我们还是要来可视化一下。用的是spplot函数,这里就不展开讲了,只贴出图(当然只是随手画的,色带啥的都没调)。

4.ArcGIS中OLS工具回归结果可视化(R语言版)

最后的最后。对应虾神文章的Python版本ArcGIS中OLS工具回归结果可视化,写个R语言版本。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#载入包
#如果没安装,请先安装,如果已安装,请注释
#install.packages(".../arcgisbinding_1.0.0.125.zip", repos = NULL, type = "win.binary")…表示arcgisbinding离线包的路径
#install.packages("car")
#install.packages("GGally")
#install.packages("ggplot2")
library(arcgisbinding)
library(car)
library(GGally)
library(ggplot2)
#设置工作路径
setwd("F:/R/demo/readdata")
#检查ArcGIS产品许可
arc.check_product()
#读取数据并将数据转换为数据框
olsdata<-arc.open("china.gdb/olstest")
olsdata
olsdataframe<-arc.select(olsdata,fields = c("gdp","Index_2000","Pop_Urban","POPU","PRODUCT","Estimated","Residual","StdResid"))
#把因变量和自变量单独分离出来并用car包里的spm函数绘图
variableframe<-olsdataframe[,c(1:5)]
spm(variableframe,diagonal="hist")

感觉似乎不是很好看,换个方式。

1
2
#利用GGally的ggpairs函数画图
ggpairs(variableframe,upper = list(continuous="cor"),lower = list(continuous="smooth"),diag = list(continuous="barDiag"))

1
2
3
4
5
#绘制标准残差的分布,用ggplot2画图
a<-ggplot(olsdataframe,aes(x=StdResid))+
geom_histogram(aes(y=..density..),binwidth = 0.5,colour="white",fill="grey")+
geom_line(stat='density',colour="#FF6666")
a

1
2
3
4
5
6
7
8
9
10
11
#绘制标准残差和观测值的散点图
opar<-par(no.readonly = T)
par(fig=c(0,0.8,0,0.8))
plot(olsdataframe$gdp,olsdataframe$StdResid,col="grey",pch=16)
par(fig=c(0,0.8,0.7,1),new=T)
hist(olsdataframe$gdp,col="grey")
par(fig=c(0.75,1,0,0.8),new=T)
hist(olsdataframe$StdResid,col="grey")

主要是为了和虾神最后的效果类似,事实上,在读取完数据框之后,纯属散点图矩阵可视化方面的内容。
最后贴出虾神的公众号和博客。

微信公众号:虾神daxialu——以推广空间分析和空间数据挖掘为己任,致力于在GIS界传递分析价值。

虾神博客原文地址:《白话空间统计二十三:回归分析番外-ArcGIS中的OLS(三)》

http://blog.csdn.net/allenlu2008/article/details/70456024

坚持原创技术分享,您的支持将鼓励我继续创作!