GDAL线面互转换（2）

     在上一个文章中介绍了线转化为面和面转化为线，其主要的实现思路就是把面中的点取出来构成线，把线中的点取出来构成面，实际上就是一个硬拷贝，无奈客户的实际需求并非如此，客户想要线转面的时候几条相交线构成面，面转线的时候相同的线去除，所以又重新对功能进行了调整。

     关于线构面，这个过程中也是挺曲折的，在GDAL的群里问了好久，大家给的答案一致是需要自己写，线构面的算法需要使用左转或者右转算法，并且在网上查了相关的资料，发现蒋波涛编著的《插件式GIS应用框架的设计与实现：基于C#和AE 9.2》一书中有ArcGIS实现GDAL的完整代码，没办法，埋头去干吧，先做了一部分前期的验证，验证的时候，公司一大牛问我OGRGeometry中的Polygonize ()方法是干啥的？仔细查阅相关资料，就是要做线转面的啊！赶紧验证，发现失败，后来又仔细分析，似乎该方法需要OGRGeometry的空间拓扑关系正确，如何保证这点儿呢？考虑用Union来把一根根线合并了拓扑关系应该没有问题吧，实际测试了一下，的确正确，心中窃喜。后来拿实际的省界线来合成省面，发现效率超级慢，跟踪代码发现，OGRGeometry合并的时候随着合并线段的增加时间也在快速的增加……问题找到，上网查找解决方案，看到两篇文章（http://www.docin.com/p-1155026547.html）（http://www.docin.com/p-1238395230.html）介绍说可以考虑使用UnionCascaded（级联求并）可以大大的加快效率，实际验证确不知道怎么使用，考虑到项目的进度要求，这个问题暂且搁置了，后续再提高吧。

     对于面构线，就相对简单了不少，先把一个个的面转化成线（硬拷贝，参考上一博文中的方法），再把转化成的线Union了，生成一个拓扑关系正确的OGRMultiPolygon，再调用Simplify方法，得到的结果即为想要的线，不过当线特别多还很复杂的时候效率也高，问题和线构面的时候一样。

     在此补充上源代码，在此留个记录，后续考虑怎么解决这个问题吧……

  1 /*
  2 * @brief ConvertPolygonToPolyline        面图层转换为线图层
  3 * @param[in] tString polylinePath        转换后线图层文件路径
  4 * @param[in] OGRLayer* pLayer            要转换的面图层文件    
  5 * @param[in] Envelope envelope            要转换数据的范围
  6 * @param[in] vector<long> vecFIDs        选中的要素ID列表
  7 * @param[in] pOGRSpatialReference        要转换的数据的空间参考(如果为空表示坐标系信息不变)
  8 * @return bool                    是否成功
  9 * @author 
 10 * @date 
 11 * @note 2015年11月04日 小八创建；
 12 */
 13 bool FeatureLayerOperator::ConvertPolygonToPolylineEx(tString polylinePath,OGRLayer* pLayer,Envelope envelope,vector<long> vecFIDs,OGRSpatialReference* pOGRSpatialReference)
 14 {
 15     // 判断
 16     if(pLayer==NULL) return false;
 17     if(true==polylinePath.empty()) return false;
 18 
 19     // 坐标系读取
 20     OGRSpatialReference* pOGRSpatialReference_Source=pLayer->GetSpatialRef();    
 21     bool isSameCoordSystem = false;
 22     if(pOGRSpatialReference == NULL)
 23     {
 24         pOGRSpatialReference=pOGRSpatialReference_Source;
 25         isSameCoordSystem=true;
 26     }
 27     else if(pOGRSpatialReference!=NULL && pOGRSpatialReference_Source!=NULL)
 28     { 
 29         isSameCoordSystem=pOGRSpatialReference_Source->IsSame(pOGRSpatialReference);
 30     }
 31 
 32     // 创建Shape文件
 33     OGRDataSource* pOGRDataSource=CreateShapeFile(polylinePath,pOGRSpatialReference,wkbLineString);
 34     if(pOGRDataSource==NULL) return false;
 35 
 36     OGRLayer* pOGRLayer=pOGRDataSource->GetLayer(0);
 37     if(pOGRLayer==NULL) return false;
 38 
 39     // 面转线再合并
 40     OGRFeature* pOGRFeature_Old;
 41     OGRGeometry* pTempGeometry=NULL;
 42     OGRGeometry* pTempGeometryUnion=NULL;
 43 
 44     // 当前选择导出
 45     if(false==vecFIDs.empty()&&vecFIDs.size()>0)
 46     {
 47         for(int i=0;i<vecFIDs.size();i++)
 48         {
 49             pOGRFeature_Old=pLayer->GetFeature(vecFIDs[i]);
 50             pTempGeometry=ConvertPolygonToPolylineGeo(pOGRFeature_Old->GetGeometryRef());
 51             pTempGeometry->assignSpatialReference(pOGRSpatialReference_Source);
 52             if(false == isSameCoordSystem)pTempGeometry->transformTo(pOGRSpatialReference);
 53 
 54             if(pTempGeometryUnion==NULL) pTempGeometryUnion=pTempGeometry;
 55             else pTempGeometryUnion=pTempGeometryUnion->Union(pTempGeometry);
 56         }
 57     }
 58     else
 59     {
 60         if(false!=envelope.isNull()&&envelope.getMaxX()!=envelope.getMinX()) 
 61         {        
 62             pLayer->SetSpatialFilterRect(envelope.getMinX(),envelope.getMinY(),envelope.getMaxX(),envelope.getMaxY());
 63         }
 64 
 65         pOGRFeature_Old=pLayer->GetNextFeature();
 66         while(NULL!= pOGRFeature_Old)
 67         {
 68             pTempGeometry=ConvertPolygonToPolylineGeo(pOGRFeature_Old->GetGeometryRef());
 69             pTempGeometry->assignSpatialReference(pOGRSpatialReference_Source);
 70             if(false == isSameCoordSystem)pTempGeometry->transformTo(pOGRSpatialReference);
 71 
 72             if(pTempGeometryUnion==NULL) pTempGeometryUnion=pTempGeometry;
 73             else pTempGeometryUnion=pTempGeometryUnion->Union(pTempGeometry);
 74 
 75             pOGRFeature_Old=pLayer->GetNextFeature();
 76         }
 77     }
 78 
 79     // 获得Simply的距离
 80     double distanceValue=0.0;
 81     if(NULL==pOGRSpatialReference) // 如果为空
 82     {
 83         OGREnvelope pTempOGREnvelope ;
 84         pTempGeometryUnion->getEnvelope(&pTempOGREnvelope);
 85         if(pTempOGREnvelope.MaxX<180) distanceValue=0.00000001;
 86         else distanceValue=0.01;
 87     }
 88     else if(true==pOGRSpatialReference->IsProjected()) // 如果是Project的
 89     {
 90         distanceValue=0.01;
 91     }
 92     else // 如果是Geo的
 93     {
 94         distanceValue=0.00000001;
 95     }
 96     pTempGeometryUnion=pTempGeometryUnion->Simplify(distanceValue);
 97 
 98     // 在ShapeFile文件中添加数据行
 99     OGRFeature* pOGRFeature_New;
100     OGRGeometry* pOGRGeometry;
101     OGRFeatureDefn* pOGRFeatureDefn=NULL;
102     pOGRFeatureDefn=pOGRLayer->GetLayerDefn();
103 
104     OGRwkbGeometryType ogrGeometryType=pTempGeometryUnion->getGeometryType();
105     ogrGeometryType=wkbFlatten(ogrGeometryType);
106 
107     if(ogrGeometryType==OGRwkbGeometryType::wkbMultiLineString)
108     {
109         OGRGeometryCollection* pOGRGeometryCollectionTarget=(OGRGeometryCollection*) pTempGeometryUnion;
110         int geometryCount=pOGRGeometryCollectionTarget->getNumGeometries();
111         for(int i=0;i<geometryCount;i++)
112         {
113             pOGRFeature_New=OGRFeature::CreateFeature(pOGRFeatureDefn);
114             pOGRGeometry=pOGRGeometryCollectionTarget->getGeometryRef(i);
115             pOGRFeature_New->SetGeometry(pOGRGeometry);
116             pOGRLayer->CreateFeature(pOGRFeature_New);
117 
118             OGRFeature::DestroyFeature(pOGRFeature_New);
119             pOGRFeature_New=NULL;
120         }
121     }
122     else if(ogrGeometryType==OGRwkbGeometryType::wkbLineString)
123     {
124         pOGRFeature_New=OGRFeature::CreateFeature(pOGRFeatureDefn);
125         pOGRGeometry=pTempGeometryUnion;
126         pOGRFeature_New->SetGeometry(pOGRGeometry);
127         pOGRLayer->CreateFeature(pOGRFeature_New);
128 
129         OGRFeature::DestroyFeature(pOGRFeature_New);
130         pOGRFeature_New=NULL;
131 
132     }
133     OGRDataSource::DestroyDataSource(pOGRDataSource);
134 
135     // 销毁pTargetGeometrys
136     OGRGeometryFactory::destroyGeometry(pTempGeometryUnion);
137     pTempGeometryUnion=NULL;
138 
139     return true;
140 }
141 
142 /*
143 * @brief ConvertPolygonToPolyline        线图层转换为面图层
144 * @param[in] tString polylinePath        转换后面图层文件路径
145 * @param[in] OGRLayer* pLayer            要转换的线图层文件    
146 * @param[in] Envelope envelope            要转换数据的范围
147 * @param[in] vector<long> vecFIDs        选中的要素ID列表
148 * @param[in] pOGRSpatialReference        要转换的数据的空间参考(如果为空表示坐标系信息不变)
149 * @return bool                    是否成功
150 * @author 
151 * @date 
152 * @note 2015年11月04日 小八创建；
153 */
154 bool FeatureLayerOperator::ConvertPolylineToPolygonEx(tString polylinePath,OGRLayer* pLayer,Envelope envelope,vector<long> vecFIDs,OGRSpatialReference* pOGRSpatialReference)
155 {
156     // 判断
157     if(pLayer==NULL) return false;
158     if(true==polylinePath.empty()) return false;
159 
160     // 坐标系读取
161     OGRSpatialReference* pOGRSpatialReference_Source=pLayer->GetSpatialRef();    
162     bool isSameCoordSystem = false;
163     if(pOGRSpatialReference == NULL)
164     {
165         pOGRSpatialReference=pOGRSpatialReference_Source;
166         isSameCoordSystem=true;
167     }
168     else if(pOGRSpatialReference!=NULL && pOGRSpatialReference_Source!=NULL)
169     { 
170         isSameCoordSystem=pOGRSpatialReference_Source->IsSame(pOGRSpatialReference);
171     }
172 
173     // 创建Shape文件
174     OGRDataSource* pOGRDataSource=CreateShapeFile(polylinePath,pOGRSpatialReference,wkbPolygon);
175     if(pOGRDataSource==NULL) return false;
176 
177     OGRLayer* pOGRLayer=pOGRDataSource->GetLayer(0);
178     if(pOGRLayer==NULL) return false;
179 
180     // 面合并
181     OGRFeature* pOGRFeature_Old;
182     OGRGeometry* pTempGeometry=NULL;
183     OGRGeometry* pTempGeometryUnion=NULL;
184 
185     // 当前选择导出
186     if(false==vecFIDs.empty()&&vecFIDs.size()>0)
187     {
188         for(int i=0;i<vecFIDs.size();i++)
189         {
190             pOGRFeature_Old=pLayer->GetFeature(vecFIDs[i]);
191             pTempGeometry=pOGRFeature_Old->GetGeometryRef();
192             if(false == isSameCoordSystem)pTempGeometry->transformTo(pOGRSpatialReference);
193 
194             if(pTempGeometryUnion==NULL) pTempGeometryUnion=pTempGeometry;
195             else pTempGeometryUnion=pTempGeometryUnion->Union(pTempGeometry);
196         }
197     }
198     else
199     {
200         if(false!=envelope.isNull()&&envelope.getMaxX()!=envelope.getMinX()) 
201         {        
202             pLayer->SetSpatialFilterRect(envelope.getMinX(),envelope.getMinY(),envelope.getMaxX(),envelope.getMaxY());
203         }
204 
205         pOGRFeature_Old=pLayer->GetNextFeature();
206         while(NULL!= pOGRFeature_Old)
207         {
208             pTempGeometry=pOGRFeature_Old->GetGeometryRef();
209             if(false == isSameCoordSystem)pTempGeometry->transformTo(pOGRSpatialReference);
210 
211             if(pTempGeometryUnion==NULL) pTempGeometryUnion=pTempGeometry;
212             else pTempGeometryUnion=pTempGeometryUnion->Union(pTempGeometry);
213 
214             pOGRFeature_Old=pLayer->GetNextFeature();
215         }
216     }
217 
218     OGRGeometry* pOGRGeometryUnion = pTempGeometryUnion->Polygonize();
219 
220     // 在ShapeFile文件中添加数据行
221     OGRFeature* pOGRFeature_New;
222     OGRGeometry* pOGRGeometry;
223     OGRFeatureDefn* pOGRFeatureDefn=NULL;
224     pOGRFeatureDefn=pOGRLayer->GetLayerDefn();
225 
226     OGRwkbGeometryType ogrGeometryType=pOGRGeometryUnion->getGeometryType();
227     ogrGeometryType=wkbFlatten(ogrGeometryType);
228 
229     if(ogrGeometryType==OGRwkbGeometryType::wkbGeometryCollection||ogrGeometryType==OGRwkbGeometryType::wkbMultiPolygon)
230     {
231         OGRGeometryCollection* pOGRGeometryCollectionTarget=(OGRGeometryCollection*) pOGRGeometryUnion;
232         int geometryCount=pOGRGeometryCollectionTarget->getNumGeometries();
233         for(int i=0;i<geometryCount;i++)
234         {
235             pOGRFeature_New=OGRFeature::CreateFeature(pOGRFeatureDefn);
236             pOGRGeometry=pOGRGeometryCollectionTarget->getGeometryRef(i);
237             pOGRFeature_New->SetGeometry(pOGRGeometry);
238             pOGRLayer->CreateFeature(pOGRFeature_New);
239 
240             OGRFeature::DestroyFeature(pOGRFeature_New);
241             pOGRFeature_New=NULL;
242         }
243     }
244     else if(ogrGeometryType==OGRwkbGeometryType::wkbPolygon)
245     {
246         pOGRFeature_New=OGRFeature::CreateFeature(pOGRFeatureDefn);
247         pOGRGeometry=pOGRGeometryUnion;
248         pOGRFeature_New->SetGeometry(pOGRGeometry);
249         pOGRLayer->CreateFeature(pOGRFeature_New);
250 
251         OGRFeature::DestroyFeature(pOGRFeature_New);
252         pOGRFeature_New=NULL;
253 
254     }
255     OGRDataSource::DestroyDataSource(pOGRDataSource);
256 
257     // 销毁pTargetGeometrys
258     OGRGeometryFactory::destroyGeometry(pTempGeometryUnion);
259     pTempGeometryUnion=NULL;
260 
261     return true;
262 }
View Code

 

posted @ 2015-11-13 16:59  松山仪龙  阅读(2320)  评论(0)  编辑  收藏  举报