GDAL1.11版本对SHP文件索引加速测试

GDAL库中对于矢量数据的读取中可以设置一些过滤器来对矢量图形进行筛选，对于Shapefile格式来说，如果数据量太大，设置这个过滤器时间慢的简直无法忍受。好在GDAL1.10版本开始支持读取Shapefile文件的空间索引文件(.sbn / .sbx)来进行加速。下面就同样的数据同样的代码来对GDAL1.9.0和GDAL1.11.0两个版本进行测试时间，比较下速度（看到结果你肯定会张大嘴巴的～～）。

首先是测试代码，功能很简单，两个shp文件，一个点文件，一个面文件。面文件很大，需要根据点文件中的点来查询到对应的面文件中的图形。在此感谢“午夜风”提供的数据进行测试。图1是使用ArcMap打开两个数据显示的效果，图2是两个数据的数据量以及要素个数。

图1 使用ArcMap打开的效果

图2 两个数据的数据量和要素个数

下面是测试代码，只贴出来关键部分的函数。

void SearchSampleDataFromSHP_liml()
{
	const char* pszPoints = "C:\Users\LiMinlu\Desktop\SHP\C5Pointnew.shp";
	const char* pszPolygs = "C:\Users\LiMinlu\Desktop\SHP\C5.shp";

	// 注册驱动以及配置项
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");
	CPLSetConfigOption("SHAPE_ENCODING","");
	OGRRegisterAll();

	//打开两个数据
	OGRSFDriver* poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");

	OGRDataSource* pPntDS = poDriver->Open(pszPoints, false);
	if (pPntDS==NULL)
	{
		cout<<"打开文件：" <<pszPoints << "失败"<<endl;
		return ;
	}

	OGRDataSource* pPlgDS = poDriver->Open(pszPolygs, false);
	if (pPlgDS==NULL)
	{
		cout<<"打开文件：" <<pszPolygs << "失败"<<endl;
		return ;
	}

	OGRFeature *pPntFeature = NULL, *pPlgFeature = NULL;

	OGRLayer* pPntLayer = pPntDS->GetLayer(0);
	pPntLayer->ResetReading();
	pPntFeature = pPntLayer->GetFeature(0);

	int nFeildCount = pPntFeature->GetFieldCount();
	int nFeatureCount = pPntLayer->GetFeatureCount();

	OGRLayer* pPlgLayer = pPlgDS->GetLayer(0);
	pPlgLayer->ResetReading();
	pPlgFeature = pPlgLayer->GetFeature(0);

	int nPlgFeildCount = pPlgFeature->GetFieldCount();

	for (int i=0;i<nFeatureCount; i++)
	{
		pPntFeature = pPntLayer->GetFeature(i);
		double dValue = pPntFeature->GetFieldAsDouble(nFeildCount-1);
		OGRPoint *pPoint = (OGRPoint *)pPntFeature->GetGeometryRef();

		//设置面图层的过滤属性
		pPlgLayer->ResetReading();
		pPlgLayer->SetSpatialFilter((OGRGeometry*)pPoint); 

		pPlgFeature = pPlgLayer->GetNextFeature();

		if(pPlgFeature == NULL)
		{
			OGRFeature::DestroyFeature(pPntFeature);
			continue;
		}

		OGRFeature::DestroyFeature(pPntFeature);
		OGRFeature::DestroyFeature(pPlgFeature);
	}

	OGRDataSource::DestroyDataSource(pPntDS);
	OGRDataSource::DestroyDataSource(pPlgDS);
}

下面是main函数以及输出时间的一个小函数。

void ShowTime()
{
	time_t t = time(0); 
	char tmp[64]; 
	strftime( tmp, sizeof(tmp), "%Y/%m/%d %X", localtime(&t) ); 
	puts( tmp ); 
}

int _tmain(int argc, _TCHAR* argv[])
{
	ShowTime();
	SearchSampleDataFromSHP_liml();
	ShowTime();

	system("pause");
	return 0;
}

首先来看看使用GDAL1.9.0版本的时间，处理时间如图3所示。（注意，以下测试时间全部使用Release版本进行测试所得）

图3 使用GDAL1.9.0所用时间

再看看GDAL1.11.0所用的时间，处理时间如图4所示。

图4 使用GDAL1.11.0所用时间

由上面两个处理时间可以看到，在GDAL1.11.0版本处理时间大幅度提高（100倍啊），所以用到了空间索引这块的同学还是将GDAL的版本更新一下吧。

我们知道shapefile文件一般必须的是3个文件，后缀名是shp、dbf和shx。如果数据有投影信息的话再加一个prj文件。这种标准的shp文件是我们常用的，使用GDAL创建的话也会生成这么几个文件。但是当用ArcMap打开的时候，会自动多出来几个文件，后缀名是sbn和sbx，可能还有个xml文件。如图2所示。这两个文件就是ArcMap自动生成的空间索引文件（ESRI spatial index files）。

按照GDAL的官方文档说明，目前GDAL库只支持读取空间索引文件，还不支持创建，所以如果要处理大数据量的shp文件，可以先用ArcMap打开让其创建好空间索引文件再用GDAL处理。此外GDAL还支持读写UMN MapServer使用的四叉树索引文件(.qix)。具体可以参考GDAL官网中的Shapefile格式页面（网址是：http://www.gdal.org/drv_shapefile.html）。