最近在使用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可以正常打开,程序也不会出现崩溃的情况,但是属性值显示仍然有点小问题,不过已经不影响使用了。