使用Python函数式编程进行空间回归分析

200

9.06No.3

WWW.GWN.O R https://www.sodocs.net/doc/611077580.html,

地理信息世界

GEOMATICS WORLD

Abstract:As the traditional programming languages do not suite for spatial data analysis,

functional programming is more advanced and more practical.This paper put forward to intro －duce Python ’s functional programming into spatial regression analysis.It analyzed and dis －cussed application prospect of script language combined with functional programming method in GIS data analysis.The experiments proved that it can effectively improve the programming

efficiency,reduce the codes,and focus the researcher ’

s attention on the problem using the Python ’s functional programming combined with the existed open source software packages.Keywords:Python;functional programming;Spatial Regression Analysis;Open Source Soft －ware;data mining

0引言

从目前涉及空间数据分析的主要教材、书籍和论文来看，对于空间分析的实现和解算，一部分研究人员采用封装了统计功能的GIS ，RS 软件模块进行操作[1]

；另一

部分人员采用诸如SPSS ，M atlab 等统计、数学计算软件实现[2]

。既有的

软件只包含有限的功能，可定制性不强；数学软件不易操作GIS 特有的数据结构，数据准备的过程较为烦琐。而在学校教学中，学生多数使用C 语言等传统的编程语言完成空间分析程序的实现，过程烦琐。对于较为灵活的探索性空间数据分析来说，需要不断尝试各种模

型进行比较分析，因此要求编程工具具备灵活性和强大数值计算能力。Python 是一种面向对象的动态语言，采用Python 进行函数式编程恰好可以解决此类问题，在编程效率和可定制性之间找到一个平衡，也便于空间分析模块整合在更复杂的程序中。

同时，Python 平台支持从数据库连接到Web 服务器的全面支持，以Python 开发的程序模块可以和系统整合参与实际应用。

1编程工具简介

Python 是一种底层基于C 语

言实现的脚本语言，诞生于1991年。作为一种脚本语言，Python 采用解释运行的方式，编写后无需进行

编译即可直接通过解释器执行，具有典型的动态语言特点，编程效率极高。Python [3]目前版本为2.6，其已从Linux 平台上的一个脚本语言发展成为跨操作系统、多种底层实现的通用语言，

并且可以实现对既有的C/C++库调用，大大丰富了Pyt-hon 的功能。

函数式编程（Functional Pro-gramming

）是一种编程范式[4]。以函数为程序基础，通过向函数传入参数，返回输出而获得计算结果，函数本身由多个函数组成，直至函数的最底层仅仅采用程序的基本运算功能实现[5]。采用函数为基础，抛弃了传统过程式语言的控制流和状态变量，对于复杂的序列计算，

使用Python 函数式编程进行空间回归分析

文章编号：1672-1586（2009）02-0077-04中图分类号：TP311.132.4文献标识码：A

摘要：针对使用传统的编程方式进行空间回归分析过于繁杂的问题，本文提出将Python 函数式编程引入空间回归分析方法，分析和探讨脚本语言+函数式编程方式在GIS 数据分

析中的推广应用前景。编程实践证明，采用Python 函数式编程，结合既有的开源软件包，可以有效提高编程效率，减少代码量，将实践人员的注意力集中在问题本身。关键词：Python ；函数式编程；空间回归分析；开源软件；数据挖掘

孙

宁，蒲英霞

（南京大学地理与海洋科学学院，江苏南京210093）

SUN Ning ，PU Ying-xia

（School of Geographic and Oceanographic Sciences ，Nanjing University ，Nanjing 210093，China ）

Spatial Regression Analysis with Python ’s Functional

Programming

孙宁（1987-），男，

江苏南京人，本科学生，南京大学地理与海洋科学学院地理信息科学系，

研究方向：

GIS 开放标准及开源软件应用，脚本语言发展及应用。基金项目：

国家自然科学

基金项目资助（40601074

）。E-mail ：classicning

@https://www.sodocs.net/doc/611077580.html,

收稿日期：2008-06-28

66

2009.0

6No.3

WWW.GWN.O R https://www.sodocs.net/doc/611077580.html, 地理信息世界

GEOMATICS WORLD

有效提高了编程效率。实现了函数式编程的语言包括Python ，Haskell ，Scheme ，Erlang 和各种数学软件专用语言（如S ，M L 等）。

Python 函数式编程支持[6]

：

1）匿名函数lambda

在Python 中，可以使用lambda 关键字定义匿名函数，并将此函数作为参数传入其他函数中。

lambda argument1，argument2，...argumentN :expression using argu-ments

2）map （）

映射函数，可以按照参数中指定的函数对序列进行映射，生成新的序列。例如：

>>>map （lambda x:x*2，range （5））

[0，2，4，6，8]3）reduce （）

reduce 操作通过指定的函数对序列进行迭代运算。例如：

>>>reduce （lambda x ，y :x ·（y +1），range （9），1）

3628804

）filter （）过滤器用户从一个序列中取出满足条件的子序列。例如：>>>filter （lambda x :x %2==0，range （10））

[0，2，4，6，8]

2空间回归分析

空间数据分析（Spatial Data

Analysis ，SDA ）是GIS 进行分析、决策的重要工具。利用空间数据分析技术，探索和挖掘现象分布的空间关系是SDA 的核心所在。SDA 的实现分为以下四个层次：①认知，对空间数据进行有效获取和科学的组织描述，利用空间数据来再现事

物本身；②解释，理解并解释空间数据的背景过程，认识事件的本质

规律；③预报，了解、掌握事件发生的规律后，运用预测模型对未来的状况做出合理推测；④宏观决策和调控，根据SDA 结果做出合理决

策，调控地理空间上发生的事件[7]

。空间数据分析通过对已存在的现象匹配最为接近的地理空间模型，建立可进行判别和预测的现象模型。空间回归分析，是将传统线性回归分析结合空间权重矩阵，从而引入空间依赖的一种分析方法。由于一部分地理现象在空间上的分

布是连续的，传统的线性回归模型针对离散的数据进行方程的拟合，没有对数据之间的关联影响进行处理，因而计算出的模型具有一定的局限性。选择适当的空间关联模型，有助于地理现象关联模式的探索，结合空间依赖关系的回归模型，本身保持了一定的连续性，因此在描述地理现象上具有更大的优势。

针对空间滞后性的特点，主要有三种空间滞后模型，划分的标准即在于空间滞后性的出现主要发生在因变量还是自变量[8]。本文以自变量的空间回归模型为例，假设

空间滞后出现在自变量X 上：

Y=X β1+WX β2+ε

式中，X 为k ×n 的自变量矩阵，W 为进行了行标准化的n ×n 空间

权重矩阵，β为k ×1的系数向量，ρ是高阶模型的空间自回归系数，ε是n ×1回归误差向量。

3空间回归分析模型解算

下面以全国各省为单元，解算各省一年内主要工业资源常量（煤炭、电力、钢铁）与工业产值关系的

空间回归模型。

3.1统一量纲

在进行回归分析计算之前，由

于各项指标的量纲不同，需要统一进行数据变换，如中心化变换、极差标准化变换、

规格化变换等。规格化变换的公式如下：

X 'ij =X ij -min{X ij }

max{X ij }-min{X ij }

X ij

表示第i 行、第j 列数据，

min{X ij }，max{X ij }分别表示第j 列最小值和最大值。def f （m ）:

return [[（float （m[i ][j ]）-min （zip

（*m ）[j ]））/（max （zip （*m ）[j ]）-min （zip （*m ）[j ]））for j in range （len （m[i ]））]for i in range （len （m ））]

采用如下标准化变换公式，可转化为均值为0，方差为1的数据：

X '

ij =（X ij -X j ）/S j

X ij 表示第i 行、

第j 列数据，X j ，S j 分别表示第j 列数据的均值和方差。def f2（m ）:

return [[

（m[i ][j ]-sum （zip （*m ）[j ]）/float （len （zip （*m ）[j ]）））/

（sum （[（k-sum （zip （*m ）[j ]）

/float

（len （zip （*m ）[j ]）））**2for k in zip （*m ）[j ]]）/len （zip （*m ）[j ]））for j in range （len （m[i ]））]for i in range （len （m ））]

3.2建立空间权重矩阵

空间权重矩阵表达了不同空

间单元之间的空间关联关系。本文所用数据为国家基础地理信息数

67

200

9.06No.3

WWW.GWN.O R https://www.sodocs.net/doc/611077580.html,

地理信息世界

GEOMATICS WORLD

据（可从国家基础地理信息网站下载），通过python 的ogr 库结合shapely 库可以计算出各省之间的邻接关系。

用ogr 读取shp 文件，必要时需将同省的多边形进行合并，把结果存储在一个以各省名为key 的dictionary （名值对）数据结构中备用。以下是读取shp 文件放入dic-tionary 中的代码。

from osgeo import ogr #从bou2_4p.shp 中读入shp =ogr.Open （'bou2_4p.shp'）layer =shp.GetLayer （0）features =[layer.GetFeature （i ）for i in range （layer.GetFeatureCount （））]

provinces ={}for f in features:

if f.GetFieldAsString （'name'）not in provinces.keys （）:

provinces[f.GetFieldAsString （'name'）]=f

else:

geo =provinces[f.GetFieldAsS-tring （'name'）].GetGeometryRef （）

geo2=f.GetGeometryRef （）provinces[f.GetFieldAsString （'name'）].SetGeometry （geo.Union （geo2））

Python 的ogr 库是C++GDAL /OGR 库的一个包装，用于实现矢量数据的存取和简单的几何操作。

然后即可对存储在provinces 字典中的各省几何数据进行两两的空间关系运算。这里采用一种名为DE-9IM 的空间位置关系判别方法[9]。根据DE-9IM 的命名规则，多边形间的邻接关系以字符串“FF2F11212”表示。通过shapely 中

几何对象的relate 方法可以比较两图形的位置关系，将以上获得的各省几何数据进行此操作即可获得邻接关系矩阵[8]。

由ogr 的Geometry 对象到shapely 的Geometry 对象的数据交换过程可以采用WKT （Well-Known Text ）作为中间介质，如：

>>>from shapely import wkt >>>heilongjiang =wkt.loads （str （f.GetGeometryRef （）））

>>>print heilongjiang.relate （neimenggu ）

FF2F11212

这样，通过比较两两几何对象的relate 方法返回值即可确定其位置关系是否满足邻接，由此建立全国各省的几何关系邻接矩阵。

在这个工业产值和资源的模型中，仅仅依靠邻接关系建立的空间权重矩阵并不一定可靠。由于现代交通运输发达，即使不是邻接的空间对象，依然可能产生相互影响。因此可以尝试建立更加切合实际情况的距离权重矩阵。考虑到计算的简便，距离权重矩阵也采用了二值化的方式，阈值为最大距离的1/8。进而采用ogr 库读取基础地理信息，进行距离的计算，建立矩阵。

3.3全局空间自相关Moran ’s I 检验

应用全局空间自相关M oran ’s

I 检验函数如下：I =n S 0

·

n

i

Σn

j =1

Σωij

（X i

-X ）

（X j

-X ）n

i

Σ

（X i

-X ）2

若采用传统的编程方式利用上述公式，需要进行循环，多次存储状态变量，编程效率较低。若采用Python 的函数式编程，

则可直接以类似公式的形式书写代码，不必纠缠于状态变量的维护，部分代码如下：I =[]

for c in range

（3）:avg =sum （[i [c ]for i in source]）/n

total =sum （[sum （map （lambda x ，y ，z:x ×y ×z ，weight[i ]，

[source[i ][c ]-avg]·n ，

[s [c ]-avg for s in source]））for i in range （n ）]）

m =sum

（map （lambda x :（x -avg ）**2，[i [c ]for i in source]））

I.append （n /s 0*total/m ）在以上的函数式编程代码中，诸如这样复杂的计算可以在一行中计算完成。代码的表达和数学公式的表达恰好是对应的，

下面给出焦炭产量、发电量和钢铁产量三个因子在不同量纲统一方式和不同空间权重矩阵组合下的Moran ’

s I 计算结果。规格化量纲、邻接矩阵模型：[0.399070，0.287646，0.410155]

规格化量纲、距离权重矩阵模型：[0.487644，0.271192，0.412432]

标准化量纲、邻接矩阵模型：[0.399309，0.287250，0.409710]

标准化量纲、距离权重矩阵模型：[0.488278，0.270333，0.412223]

检测结果显示第四组模型的空间自相关程度最高。因此采用第四组模型进行后续的空间回归分析。

3.4空间滞后模型

本文选择标准化量纲距离权

重矩阵模型进行空间滞后模型拟合，同时采用Python 函数式编程方

68

2009.0

6No.3

WWW.GWN.O R https://www.sodocs.net/doc/611077580.html,

地理信息世界

GEOMATICS WORLD

法，结合既有的矩阵运算类库，实现系数向量的求解。

对于本模型三列的情况，可以使用以下方式实现每一行的空间权重赋值。for i in range （n ）:

wx1.append

（sum （map （lambda x ，y:x ·y ，x 1，

weight[i ]）））w x 2.append （sum （map （lambda x ，y:x ·y ，x 2，

weight[i ]）））w x 3.append （sum （map （lambda x ，y :x ·y ，x 3，

weight[i ]）））wy.append （sum （map （lambda x ，y :x ·y ，y ，

weight[i ]）））对于最小二乘法的计算模型：

β=（X T X ）-1X T

Y

那么对于自变量的空间滞后模型，X ，Y 矩阵的形式为：

X ，Y =zip （*[x 1，x 2，x 3，

w x 1，w x 2，w x 3，const]），zip （*[y ]）

这样，结合一个矩阵运算的函数库即可解算出β矩阵，得自变量的空间回归模型：

y =X

·-0.2026

0.83030.164!"""""""""#

$

%%%%%%%%%&

7

+WX ·-0.52230.08010.449!

"""""""""#

$%%%%%%%%%&

2

+0.0094

4运算脚本性能分析

对于函数式编程的性能，可以

通过Python 自带的Profile 进行运行时间的统计。

以计算M oran ’s I 的脚本为例，我们可以比较函数式编程和采用传统方法编程的效率。

在执行Python 程序时，通过python m cProfile ***.py 的方式运行，调用cProfile 模块进行Profile 统计。

对于函数式编程，Profile 结果

如下：

3440function calls （3409prim-itive calls ）in 0.019CPU seconds

对于传统的编程方式，采用循环遍历、状态变量存储，Profile 结果如下：

461function calls （430primitive calls ）in 0.009CPU seconds

函数式编程执行了更多的函数调用，在效率上有一定损失。

5结语

通过以上示例，Python 函数式编程的优缺点已经非常清晰。在较为复杂的空间分析运算中，面对大量的公式和数据，采用传统的语言编程过程过于复杂，而纯粹通过数学软件来实现需要完成一定工作

量的数据准备，数学功能相对独立，不易与既有系统衔接。Python 作为ArcGIS 9.0的使用脚本语言，在GIS 上已经得到了一定程度的推广，而其开源GIS 类库直接继承自比较成熟C++系，GDAL/OGR ，Shapely ，RTree 等都处在一个良好的发展态势。此外，除GIS 相关功能以外，Python 完全可以胜任从Web 开发到跨平台桌面开发、游戏开发等各个方面。随着脚本语言Ruby 在Web 开发上暴发式的巨大成功和Linux 操作系统的进一步推广，包括Python 在内的其他一些语言也得到了更大的关注，社区发展迅速，未来必然将会出现更多的Python 项目。

对于长期遭受诟病的解释型脚本语言的效率，尽管Python 的运行效率低于C 语言，但是随着新版本的开发，Python 的性能不断提升。另一方面，由于可以和编译后的C 模块共同使用，可将性能要求较高的部分以C 语言

实现。

Python 的最大优势在于极高的编程效率，如SDA 这样需要不断进行探索计算，进行模型的变换、参数的修改而获得最优解的研究方法，对于编程工具的敏捷性要求则更高，而这正是脚本语言的优势所在。作为一门结合了函数式编程的脚本语言，Python 的优势更加明显。因此，Python 语言在GIS 研究中的地位也必将随着人们对它的了解加深而显得越发重要。

参考文献

[1]周兰霞.ArcGIS9.0在苯丙酮尿症空间

分析中的应用[J].中国卫生统计，2007，24（6）:604.

[2]刘承良.武汉都市圈经济社会要素流

的空间分析[J].人文地理，2007，22（6）:30.

[3]Python Official Documentations[EB/OL].

2008-04-16.https://www.sodocs.net/doc/611077580.html,/a-bout/.

[4]Wikipedia:Functional Programming.

[EB/OL].2008-04-16.https://www.sodocs.net/doc/611077580.html,/wiki/Functional_programming [5]John Hughes.Why Functional Program-ming M atters [EB/OL].2008-04-18.http://www.cschalmers.se/~rjmb/papers/vhy-fp.html.

[6]M ark Lutz.Learning Python 3rd Edition

[M ].O ’Reilly ，2007:350-353.[7]朱欣娟.基于GIS 的空间分析及其发

展研究[J].计算机工程与应用2002，18（1）:1.

[8]L.Anselin.SPATIAL ECONOM ETRICS

[M ].Springer ，1988:5-11.

[9]OpenGIS Implementation Specification

for Geographic information -Simple fea-ture access -Part 1:Common architec-ture [EB/OL].2008-03-20.https://www.sodocs.net/doc/611077580.html,/standasds/sfa.

69