GDAL不支持创建PCIDSK的面状矢量格式

最近在使用GDAL创建PCIDSK格式的矢量数据，发现创建点和线的矢量数据都没问题，创建面状的只有属性表没有图形。在GDAL官网说明也写的是支持的，地址为：http://www.gdal.org/frmt_pcidsk.html。

实在没办法，翻看GDAL源码才发现，SetFeature的时候，只写了wkbPoint和wkbLineString类型，其他的加了句Debug代码如下：

CPLDebug( "PCIDSK", "Unsupported geometry type in SetFeature(): %s",
poGeometry->getGeometryName() );

这也太……没办法只好自己研究研究补齐了。

通过查看GDAL读取发现，对于面状的矢量，PCIDSK在矢量段里面的属性表里面多存了一个字段，叫RingStart，类型为IntList，用来存储多边形中各个环的起始点号。对于PCIDSK的矢量数据，所有的点都是存储在一个大的数组里面的，对于多边形而言，可能有多个环组成，这些环里面所有的点全部都存在一个数组中，如何来区分每个环的顶点坐标，就需要通过RingStart里面的值来进行分割。知道了存储的原理，那么就按照这个原理将写多边形的部分补上就可以了。修改后的代码如下：

else if( wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon )
{
OGRPolygon *poPoly = (OGRPolygon *) poGeometry;
int nRingSize = poPoly->getNumInteriorRings();

std::vector<PCIDSK::int32> anRingStart;
anRingStart.resize(nRingSize+1);

OGRLinearRing *poRing = NULL;
poRing = poPoly->getExteriorRing();
anRingStart[0] = poRing->getNumPoints();

aoVertices.resize(poRing->getNumPoints());
for(int i = 0; i < aoVertices.size(); i++ )
{
aoVertices[i].x = poRing->getX(i);
aoVertices[i].y = poRing->getY(i);
aoVertices[i].z = poRing->getZ(i);
}

for (int iRing=0; iRing < nRingSize; iRing++)
{
int nCurrentStart = aoVertices.size();
poRing = poPoly->getInteriorRing(iRing);
anRingStart[iRing+1] = nCurrentStart + poRing->getNumPoints();
aoVertices.resize(nCurrentStart + poRing->getNumPoints());

for(int i = nCurrentStart; i < aoVertices.size(); i++ )
{
aoVertices[i].x = poRing->getX(i-nCurrentStart);
aoVertices[i].y = poRing->getY(i-nCurrentStart);
aoVertices[i].z = poRing->getZ(i-nCurrentStart);
}
}

if(iRingStartField == -1)
{
iRingStartField = poVecSeg->GetFieldCount();
//poVecSeg->AddField( "RingStart", PCIDSK::FieldTypeCountedInt, "", "" );
OGRFieldDefn oField( "RingStart", OFTIntegerList );
//oField.SetWidth(100);
CreateField( &oField );
}

std::vector<PCIDSK::ShapeField> aoShapeFields;
poVecSeg->GetFields(id, aoShapeFields);

aoShapeFields[iRingStartField].SetValue(anRingStart);
poVecSeg->SetFields( id, aoShapeFields );
}

修改完之后，重新编译GDAL就可以了。

=======================修改于2015年1月9日========================

通过上面的代码是可以生成一个面状的pix文件，使用GDAL打开也没问题，但是使用Focus或者MosaicTool打开的时候会造成这两个程序崩溃，同时属性值显示会有一定点问题。今天再仔细查看了下读取部分的代码，发现对于没有内环的多边形，也就是说一个Feature里面只有一个多边形的时候，不需要写RingStart这个属性值，只有含油内环的时候，也就是多边形中有内环的时候才需要，所以将上面的代码修改为下面的代码：

else if( wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon )
{
OGRPolygon *poPoly = (OGRPolygon *) poGeometry;
OGRLinearRing *poRing = NULL;
poRing = poPoly->getExteriorRing();

aoVertices.resize(poRing->getNumPoints());
for(int i = 0; i < aoVertices.size(); i++ )
{
aoVertices[i].x = poRing->getX(i);
aoVertices[i].y = poRing->getY(i);
aoVertices[i].z = poRing->getZ(i);
}

int nRingSize = poPoly->getNumInteriorRings();
if(nRingSize > 0 )
{
std::vector<PCIDSK::int32> anRingStart;
anRingStart.resize(nRingSize+1);
anRingStart[0] = poRing->getNumPoints();

for (int iRing=0; iRing < nRingSize; iRing++)
{
int nCurrentStart = aoVertices.size();
poRing = poPoly->getInteriorRing(iRing);
anRingStart[iRing+1] = nCurrentStart + poRing->getNumPoints();
aoVertices.resize(nCurrentStart + poRing->getNumPoints());

for(int i = nCurrentStart; i < aoVertices.size(); i++ )
{
aoVertices[i].x = poRing->getX(i-nCurrentStart);
aoVertices[i].y = poRing->getY(i-nCurrentStart);
aoVertices[i].z = poRing->getZ(i-nCurrentStart);
}
}

if(iRingStartField == -1)
{
iRingStartField = poVecSeg->GetFieldCount();
OGRFieldDefn oField( "RingStart", OFTIntegerList );
CreateField( &oField );
}

std::vector<PCIDSK::ShapeField> aoShapeFields;
poVecSeg->GetFields(id, aoShapeFields);
aoShapeFields[iRingStartField].SetValue(anRingStart);
poVecSeg->SetFields( id, aoShapeFields );
}
}

通过测试发现，这下生成的使用Foucs和MosaicTool可以正常打开，程序也不会出现崩溃的情况，但是属性值显示仍然有点小问题，不过已经不影响使用了。