基于WebService的样方分析设计与实现

样方分析法的概念和步骤相对来说比较简单和直观。样方分析一般可分为以下四个步骤：

（1）将研究区域划分为规则的正方形网格区域；

（2）统计落入每个网格中点的数量。由于点在空间上分布的疏密性，有的网格中点的数量多，有的网格中点的数量少，有的网格中点的数量甚至为零；

（3）统计出包含不同数量点的网格数量的频率分布。

（4）将观测得到的频率分布和已知的频率分布或理论上的随机分布（如泊松分布）作比较，判断点模式的类型。

我们采用t检验方法，这部分内容参见“ArcView GIS？与ArcGIS？地理信息统计分析_（美）戴维·W·S·黄，杰·李著”一书。

既然样方分析是对研究区域内的点进行分析，那么样方分析功能的数据输入自然包括点数据和分析范围。对于点数据的输入则很简单，只需给定点shapefile数据在服务器上的绝对路径即可。但是，分析范围的输入方式就不止一种了。比如用户可以输入分析范围的左下角和右上角坐标，也可以指定一个shapefile，将其属性Extent作为分析范围，或者直接利用点shapefile的Extent属性作为分析范围等等。

?考虑到输入多样性，本文设计了三个WebMethod：

1）样方分析WebMethod1

输入：InputPointFile

????? minX,minY,maxX,maxY

2）样方分析WebMethod2

输入：InputPointFile

InputSHPFile

3）样方分析WebMethod3

输入：InputPointFile

其中，InputPointFile表示的点数据在服务器上的绝对路径；样方分析WebMethod1中minX,minY和maxX,maxY分别表示分析范围的左下角和右上角坐标；样方分析WebMethod2中InputSHPFile表示作为分析范围的SHP数据在服务器上的绝对路径；那么，样方分析WebMethod3则是利用点shapefile的Extent属性作为分析范围。

样方分析输入参数已经确定，接下来确定输出参数。三个WebMethod虽然输入参数不一样，但是输出参数是一样的。

既然是样方分析，那么输出结果中肯定要包含最终的结论，即点的分布类

型。同时样方分析过程中会产生一些文件，比如绘制的格网shapefile文件和存储统计信息的XML文件，然后为了便于下载对文件进行压缩产生的文件等。所以这些文件存放路径也需要输出，便于服务调用者后续操作。

另外，服务在运行过程中可能出现错误，导致服务非正常结束，那么输出参数中有必要包含错误信息和服务是否运行成功的信息。由于输出参数过多，本位在输出设计时，就设计了一个类，类中包含所要输出的参数，详细输出参数如下：

Class ResultMessage

{

???????? stringJobId;// 30位随机唯一编码

???????? stringErrorMessage;

???????? stringOutputStatisticPath;

???????? stringOutputGridShpPath;

???????? stringCompressFilePath;

???????? doubleValueT;

???????? intAllPnum;

???????? intAllQnum;

???????? boolIsSuccessful;

???????? stringDistributePattern;

???????? int[]Statistic;

}

这里有必要对参数进行一下说明：

JobId: 随机生成的30个字符，用于文件夹的命名

????????????? ErrorMessage:错误信息

????????????? OutputStatisticPath:包含统计信息的XML文件存放路径

????????????? OutputGridShpPath:格网存放路径

????????????? CompressFilePath:压缩后的文件存放路径

????????????? ValueT:t检验中所计算的t值

????????????? AllPnum:待分析的点总数

????????????? AllQnum:样方总数

????????????? IsSuccessful:样方分析是否成功

????????????? DistributePattern:所计算的点的分布模式

????????????? []Statistic:统计信息

本文前面已经提及，Web Service会生成XML文件。那么XML的文件结构是怎样的呢？因为XML文件记录的是统计信息，即点数和含有此点数的样方数，文件内容很简单，所以对应文件结构很简单。本文设计的XML文件结构如下：

<Statistics>

<PQuadratNum>

<Pnumber> number </Pnumber>

<Qnumber> number </Qnumber>

</PQuadratNum>

<PQuadratNum>

<Pnumber> number </Pnumber>

<Qnumber> number </Qnumber>

</PQuadratNum>

......

</Statistics>

功能设计

本文将样方分析的流程分为几个步骤，一个步骤对应一个功能。一个功能一般由一到两个函数完成。大体分为以下几个功能：

1）获取点数AllPnum

2）计算样方数AllQnum(row*col)

3）画格网Grid.shp

4）统计每个样方中的点数

5）统计含0、1、2…个点的样方数Statistic

6）计算t值判定分布模式

需要说明的是，计算样方数的时候，我们不是直接计算样方总数，而是计

算行和列上样方数。计算方法是，用分析区域边长除以最优样方边长，然后分别得到分析区域行和列上格网数，两者相乘就是样方总数。为什么不直接计算样方总数呢？因为直接计算样方总数没有意义，比如我们求得样方总数为48，那么绘制格网时，是绘制4*12的格网还是6*8的格网呢？无法根据样方总数来绘制格网。

整个功能流程图如下图所示：

GDAL是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。本文主要用到GDAL的一个分支OGR。OGR与GDAL类似，只不过它是对矢量数据的支持。

OGR所提供的功能使用起来跟使用AE类似。比如获取一个图层的要素数目可使用GetFeatureCount()方法，获取图层范围可调用其GetExtent()方法，通过空间过滤筛选要素可使用SetSpatialFilterRect()方法，等等。有了AE开发的经验使用GDAL起来很简单。但是GDAL有个缺点就是它对中文的支持很不好，所以在输入路径参数时最好不要包含中文，否则，程序多半因异常退出。

本节接下来重点讲述各个功能的实现方法。

1）获取点总数AllPnum

根据输入的点路径调用OGR读取shapefile文件的函数read()方法，然后返

回类DataSource类型的值。紧接着，调用类DataSource的GetLayerByIndex()方法获取一个图层。最后调用图层的GetFeatureCount()方法获取点总数AllPnum。

2）计算样方数AllQnum(row*col)

为了计算样方尺寸，需要知道分析区域面积。而分析区域的来源可能是坐

标，也可能是shapefile文件。如果是坐标那就很好办，直接利用坐标来计算分析区域面积，如果是shapefile文件，那么需要读取这个文件以获取其Extent属性，然后计算分析区域面积A。

?????? 获得了分析区域面积后，就可以根据公式计算最优样方尺寸L了。接下来由分析区域边长除以L得到分析区域行和列方向上样方数row和col。那么，样方数AllQnum则为row*col。这里需要注意的是，由于分析区域边长除以样方边长L不一定是整数，如果为浮点数，那么为了格网能覆盖整个研究区域需要取大于等于row和col的最小整数。

3）绘制格网Grid.shp

因为这个功能会产生文件，且要存储在服务器上。那么生成shapefile文件

之前需要判定存储位置是否存在同名的shapefile文件，如果存在同名文件则需将其删除。

可是这里就存在一个问题，shapefile文件包含dbf、shp、shx、prj等格式的多个同名文件。对shapefile进行删除的时候岂不是要对文件夹下所有文件进行判断然后决定是否删除？如果文件夹下文件很多，那么这个过程就会耗费很多时间，使得程序效率低下。幸好，GDAL帮我们解决了这样一个问题。GDAL提供了删除数据源的方法—DeleteDataSource()。只需指定shapefile文件中shp格式的文件绝对路径，即可删除整个shapefile文件。

?????? 紧接着，利用CreateLayer()方法创建图层。创建图层就需要考虑设置空间参考的问题。如果不考虑格网与点数据在一起显示的问题，我们完全不必在意格网空间参考的问题，因为其对于整个样方分析没有影响。但是，为了使格网和点叠加在一起显示能够一起显示，本文将格网的空间参考设置为与点一样的空间参考。

?????? 以上准备工作做完后，就可以开始创建正方形要素了，即样方。这部分主体程序是两个嵌套的for循环。循环里先创建一个要素，然后创建一个Geometry对象，并创建一个由5对坐标表示的正方形。Geometry对象就是基于这个正方形来创建的。这里为什么是5对坐标而不是4对呢？因为5对坐标中第一对和第五对是同一个点，表示多边形闭合了。这一点很重要。最后将要素添加到图层中，照此循环。

4）统计每个样方中的点数

统计样方中的点数可通过空间过滤的方法来筛选符合条件的点来进行。这里主要用到图层的SetSpatialFilterRect()方法，该方法中的参数就是矩形的左下角和右上角坐标。筛选出符合条件的点后，再调用图层的GetFeatureCount()方法获得符合条件的点数。如此循环。

但是，有一点需要注意。每次循环前需要对图层进行重置，否则，后续的筛选都是在前一次筛选得出的点的基础上进行筛选，即后续的筛选得出的点数全为0。

5）统计含0、1、2…个点的样方数Statistic

这一步是个比较耗时的步骤。本文采取的办法是现将步骤4中的统计结果

进行排序，排序之后包含相同点数的样方则排列在一起。紧接着，获取队列第一个值value,然后将其从队列中删除，统计量k自增一。判断value是否等于队列第一个值（此时，队列为新队列），若相等则重复上面过程，若不相等则含有相同点数的样方数为k，然后退出内循环，进行下一个外循环。最后将统计信息写到XML文件中。

6）计算t值判定分布模式

步骤5中已经统计含有0、1、2…个点对应的样方数，那么接着就应该计算

其均值和方差。均值和方差只是中间量，我们主要是计算t值。计算t值后将其与1.96进行比较。若t绝对值小于1.96则可判定点属于随机分布，若t绝对值大于1.96，这时候我们需要判断t值得符号了，若t为正，则点属于聚集分布，反之，属于分散分布（也叫平均分布）。

===========================================================================================

这里，我需要说明一下，由于自己理解有误，自己的webService代码有一个地方不是很正确。就是分析范围这个地方，有意义的范围应该是行政边界指定的范围，那么格网中与行政边界没有交集的样方就不应该参与统计（为了方便，还是要画出来），所以这一点希望读者注意。

===========================================================================================

上面格网就是程序绘制的格网，样方分析WEB服务代码， 请点击这里。

如果有疑问，欢迎大家给我留言。

调用样方分析webservice的前台界面：