src/HYDROData/HYDROData_ShapeFile.cxx

   1 // Copyright (C) 2014-2015  EDF-R&D
   2 // This library is free software; you can redistribute it and/or
   3 // modify it under the terms of the GNU Lesser General Public
   4 // License as published by the Free Software Foundation; either
   5 // version 2.1 of the License, or (at your option) any later version.
   6 //
   7 // This library is distributed in the hope that it will be useful,
   8 // but WITHOUT ANY WARRANTY; without even the implied warranty of
   9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  10 // Lesser General Public License for more details.
  11 //
  12 // You should have received a copy of the GNU Lesser General Public
  13 // License along with this library; if not, write to the Free Software
  14 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  15 //
  16 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  17 //
  18
  19 #include <HYDROData_ShapeFile.h>
  20 #include <HYDROData_PolylineXY.h>
  21 #include <HYDROData_Polyline3D.h>
  22 #include <HYDROData_Bathymetry.h>
  23 #include <HYDROData_Profile.h>
  24 #include <HYDROData_Iterator.h>
  25 #include <HYDROData_LandCoverMap.h>
  26
  27 #include <QFile>
  28 #include <QFileInfo>
  29 #include <TopoDS.hxx>
  30 #include <TopExp_Explorer.hxx>
  31 #include <TopoDS_Wire.hxx>
  32 #include <TopoDS_Vertex.hxx>
  33 #include <TopoDS_Edge.hxx>
  34 #include <TopoDS_Face.hxx>
  35 #include <BRep_Tool.hxx>
  36 #include <BRepTools.hxx>
  37 #include <Precision.hxx>
  38 #include <Geom_Curve.hxx>
  39 #include <Geom_Line.hxx>
  40 #include <Geom_TrimmedCurve.hxx>
  41 #include <Geom_TrimmedCurve.hxx>
  42 #include <BRepBuilderAPI_MakeEdge.hxx>
  43 #include <BRepBuilderAPI_MakeWire.hxx>
  44 #include <BRepBuilderAPI_MakeFace.hxx>
  45 #include <gp_Pln.hxx>
  46 #include <BRepLib.hxx>
  47 #include <ShapeFix_Shape.hxx>
  48 #include <TopTools_SequenceOfShape.hxx>
  49 #include <QColor>
  50 #include <BRepTopAdaptor_FClass2d.hxx>
  51 #include <TopExp.hxx>
  52 #include <OSD_Timer.hxx>
  53 #include <BRepLib_MakeVertex.hxx>
  54 #include <NCollection_List.hxx>
  55 #include <GC_MakeSegment.hxx>
  56 #include <BRep_Builder.hxx>
  57 #include <BRepAdaptor_Curve.hxx>
  58 #include <GCPnts_QuasiUniformDeflection.hxx>
  59 #include <TopLoc_Location.hxx>
  60 #include <Geom_Plane.hxx>
  61 #include <NCollection_Array1.hxx>
  62 #include <BRepBndLib.hxx>
  63 #include <Bnd_Box.hxx>
  64
  65 #ifdef WIN32
  66   #pragma warning( disable: 4996 )
  67 #endif
  68
  69 //SHP->TFaces (Import)
  70 #define OSD_TIMER
  71
  72 HYDROData_ShapeFile::HYDROData_ShapeFile() : myHSHP(NULL)
  73 {
  74 }
  75
  76 HYDROData_ShapeFile::~HYDROData_ShapeFile()
  77 {
  78   Free();
  79 }
  80
  81 void HYDROData_ShapeFile::Export(Handle(HYDROData_Document) theDocument,
  82                                  const QString& aFileName,
  83                                  NCollection_Sequence<Handle(HYDROData_PolylineXY)> aPolyXYSeq,
  84                                  NCollection_Sequence<Handle(HYDROData_Polyline3D)> aPoly3DSeq,
  85                                  QStringList& aNonExpList)
  86 {
  87   SHPHandle hSHPHandle;
  88   if (!aPolyXYSeq.IsEmpty() && aPoly3DSeq.IsEmpty())
  89   {
  90     hSHPHandle = SHPCreate( aFileName.toLatin1().data(), SHPT_ARC );
  91     for (int i = 1; i <= aPolyXYSeq.Size(); i++)
  92       if (WriteObjectPolyXY(theDocument, hSHPHandle, aPolyXYSeq(i)) != 1)
  93         aNonExpList.append(aPolyXYSeq(i)->GetName());
  94   }
  95   else if (aPolyXYSeq.IsEmpty() && !aPoly3DSeq.IsEmpty())
  96   {
  97     hSHPHandle = SHPCreate( aFileName.toLatin1().data(), SHPT_ARCZ );
  98     for (int i = 1; i <= aPoly3DSeq.Size(); i++)
  99       if (WriteObjectPoly3D(theDocument, hSHPHandle, aPoly3DSeq(i)) != 1)
 100         aNonExpList.append(aPoly3DSeq(i)->GetName());
 101   }
 102   if (hSHPHandle->nRecords > 0)
 103     SHPClose( hSHPHandle );
 104   else
 105   {
 106     SHPClose( hSHPHandle );
 107     QString aFN = aFileName.simplified();
 108     remove (aFN.toStdString().c_str());
 109     remove ((aFN.simplified().replace( aFN.simplified().size() - 4, 4, ".shx")).toStdString().c_str());
 110   }
 111 }
 112
 113 void HYDROData_ShapeFile::Export(Handle(HYDROData_Document) theDocument,
 114                                  const QString& aFileName,
 115                                  const Handle(HYDROData_LandCoverMap)& aLCM,
 116                                  QStringList& aNonExpList,
 117                                  bool bCheckLinear, bool bUseDiscr, double theDefl)
 118 {
 119   if (bCheckLinear && !aLCM->CheckLinear())
 120     return;
 121   //
 122   SHPHandle hSHPHandle = NULL;
 123   if ( !aLCM.IsNull() && !aLCM->IsEmpty())
 124   {
 125     hSHPHandle = SHPCreate( aFileName.toLatin1().data(), SHPT_POLYGON );
 126     HYDROData_LandCoverMap::Explorer It( aLCM );
 127     for( ; It.More(); It.Next())
 128     {
 129       TopoDS_Face aFace = It.Face();
 130       if (WriteObjectPolygon(theDocument, hSHPHandle, aFace, bUseDiscr, theDefl) != 1)
 131         aNonExpList.append(aLCM->GetName() + "_" +  QString::number(It.Index()));
 132     }
 133   }
 134   if (hSHPHandle && hSHPHandle->nRecords > 0)
 135     SHPClose( hSHPHandle );
 136   else
 137   {
 138     SHPClose( hSHPHandle );
 139     QString aFN = aFileName.simplified();
 140     remove (aFN.toStdString().c_str());
 141     remove ((aFN.simplified().replace( aFN.simplified().size() - 4, 4, ".shx")).toStdString().c_str());
 142   }
 143 }
 144
 145
 146 int HYDROData_ShapeFile::WriteObjectPolyXY(Handle(HYDROData_Document) theDocument,
 147                                            SHPHandle theShpHandle,
 148                                            Handle(HYDROData_PolylineXY) thePoly )
 149 {
 150   SHPObject     *aSHPObj;
 151   std::vector<double> x, y;
 152   std::vector<int> anPartStart;
 153
 154   for (int i = 0; i < thePoly->NbSections(); i++)
 155     if (thePoly->GetSectionType(i) == HYDROData_IPolyline::SECTION_SPLINE)
 156       return -1;
 157
 158   for (int i = 0; i < thePoly->NbSections(); i++)
 159   {
 160     anPartStart.push_back(x.size());
 161     HYDROData_PolylineXY::PointsList aPointList = thePoly->GetPoints(i);
 162     for (int j = 1; j <= aPointList.Size(); j++)
 163     {
 164       gp_XY P = aPointList(j);
 165       theDocument->Transform(P, false);
 166       x.push_back( P.X());
 167       y.push_back( P.Y());
 168     }
 169     if (thePoly->IsClosedSection(i))
 170     {
 171       gp_XY P = aPointList(1);
 172       theDocument->Transform(P, false);
 173       x.push_back( P.X());
 174       y.push_back( P.Y());
 175     }
 176   }
 177
 178   aSHPObj = SHPCreateObject( SHPT_ARC, -1, thePoly->NbSections(), &anPartStart[0], NULL, x.size(), &x[0], &y[0], NULL, NULL );
 179   SHPWriteObject( theShpHandle, -1, aSHPObj );
 180   SHPDestroyObject( aSHPObj );
 181   return 1;
 182 }
 183
 184 int HYDROData_ShapeFile::WriteObjectPoly3D(Handle(HYDROData_Document) theDocument,
 185                                            SHPHandle theShpHandle,
 186                                            Handle(HYDROData_Polyline3D) thePoly )
 187 {
 188   SHPObject     *aSHPObj;
 189   std::vector<double> x, y, z;
 190   std::vector<int> anPartStart;
 191
 192   for (int i = 0; i < thePoly->GetPolylineXY()->NbSections(); i++)
 193     if (thePoly->GetPolylineXY()->GetSectionType(i) == HYDROData_IPolyline::SECTION_SPLINE)
 194       return -1;
 195
 196   for (int i = 0; i < thePoly->GetPolylineXY()->NbSections(); i++)
 197   {
 198     anPartStart.push_back(x.size());
 199     HYDROData_PolylineXY::PointsList aPointList = thePoly->GetPolylineXY()->GetPoints(i);
 200     for (int j = 1; j <= aPointList.Size(); j++)
 201     {
 202       gp_XY P = aPointList(j);
 203       theDocument->Transform(P, false);
 204       x.push_back( P.X());
 205       y.push_back( P.Y());
 206       z.push_back(thePoly->GetAltitudeObject()->GetAltitudeForPoint(gp_XY (aPointList(j).X(), aPointList(j).Y())));
 207     }
 208     if ( thePoly->GetPolylineXY()->IsClosedSection(i))
 209     {
 210       gp_XY P = aPointList(1);
 211       theDocument->Transform(P, false);
 212       x.push_back( P.X());
 213       y.push_back( P.Y());
 214       z.push_back(thePoly->GetAltitudeObject()->GetAltitudeForPoint(gp_XY (aPointList(1).X(), aPointList(1).Y())));
 215
 216     }
 217   }
 218
 219   aSHPObj = SHPCreateObject( SHPT_ARCZ, -1, thePoly->GetPolylineXY()->NbSections(), &anPartStart[0], NULL, x.size(), &x[0], &y[0], &z[0], NULL );
 220   SHPWriteObject( theShpHandle, -1, aSHPObj );
 221   SHPDestroyObject( aSHPObj );
 222   return 1;
 223 }
 224
 225 int HYDROData_ShapeFile::WriteObjectPolygon(Handle(HYDROData_Document) theDocument,
 226                                             SHPHandle theShpHandle,
 227                                             const TopoDS_Shape& theInputShape,
 228                                             bool bUseDiscr, double theDefl)
 229 {
 230   if (theInputShape.IsNull())
 231     return 0;
 232
 233   if (theInputShape.ShapeType() == TopAbs_FACE)
 234   {
 235     ProcessFace(theDocument, TopoDS::Face(theInputShape), theShpHandle, bUseDiscr, theDefl);
 236   }
 237   else if (theInputShape.ShapeType() == TopAbs_COMPOUND)
 238   {
 239     TopExp_Explorer Ex(theInputShape, TopAbs_FACE);
 240     for (; Ex.More(); Ex.Next())
 241     {
 242       TopoDS_Face aF = TopoDS::Face(Ex.Current());
 243       if (aF.IsNull())
 244         continue;
 245       ProcessFace(theDocument, aF, theShpHandle, bUseDiscr, theDefl);
 246     }
 247   }
 248   else
 249     return 0;
 250
 251   return 1;
 252
 253 }
 254
 255 void HYDROData_ShapeFile::ProcessFace(Handle(HYDROData_Document) theDocument,
 256                                       const TopoDS_Face& theFace,
 257                                       SHPHandle theShpHandle,
 258                                       bool bUseDiscr, double theDefl )
 259 {
 260   if (theFace.ShapeType() != TopAbs_FACE)
 261      return;
 262   SHPObject     *aSHPObj;
 263   std::vector<double> x, y;
 264   std::vector<int> anPartStart;
 265   TopoDS_Wire OuterW = BRepTools::OuterWire(theFace);
 266   NCollection_Sequence<TopoDS_Wire> aWires;
 267
 268   //write an outer wire first
 269   aWires.Append(OuterW);
 270   TopExp_Explorer Ex(theFace, TopAbs_WIRE);
 271   for (; Ex.More(); Ex.Next())
 272   {
 273     TopoDS_Wire aW = TopoDS::Wire(Ex.Current());
 274     if (aW.IsEqual(OuterW))
 275       continue;
 276     aWires.Append(aW);
 277   }
 278
 279   int NbWires = 0;
 280   for (int k = 1; k <= aWires.Length(); k++)
 281   {
 282     TopoDS_Wire aW = aWires(k);
 283     if (aW.IsNull())
 284       continue;
 285     NbWires++;
 286     if (aW.Orientation() == TopAbs_INTERNAL)
 287       //cant write internal wires/edges
 288       continue;
 289     // Try to reorder edges
 290     Handle(ShapeFix_Wire) aSFW = new ShapeFix_Wire( aW, theFace, Precision::Confusion() );
 291     aSFW->ModifyTopologyMode() = Standard_False;
 292     aSFW->ModifyGeometryMode() = Standard_False;
 293     aSFW->FixReorder();
 294     Handle(ShapeExtend_WireData) aSEWD = aSFW->WireData();
 295     Standard_Integer nbE = aSEWD->NbEdges();
 296     if (nbE == 0)
 297       continue;
 298     //
 299     anPartStart.push_back(x.size());
 300     NCollection_Sequence<gp_Pnt2d> aPnts;
 301     for (Standard_Integer i = 1; i <= nbE; i++)
 302     {
 303       TopoDS_Edge E = aSEWD->Edge(i);
 304       if (!bUseDiscr)
 305       {
 306         TopoDS_Vertex aV = TopExp::LastVertex(E, 1);
 307         if (aV.IsNull())
 308           continue;
 309         gp_Pnt P = BRep_Tool::Pnt(aV);
 310         aPnts.Append(gp_Pnt2d(P.X(), P.Y()));
 311       }
 312       else
 313       {
 314         BRepAdaptor_Curve Cur( E );
 315         GCPnts_QuasiUniformDeflection Discr( Cur, theDefl );
 316         if( !Discr.IsDone() )
 317           continue; //skip edge?
 318         double NewDefl = theDefl/2.0;
 319         while (Discr.NbPoints() < 2)
 320         {
 321           Discr.Initialize(Cur, NewDefl);
 322           NewDefl = NewDefl/2.0;
 323         }
 324         //
 325         if (E.Orientation() == TopAbs_FORWARD)
 326         {
 327           for( int i = 1; i <= Discr.NbPoints(); i++ )
 328           {
 329             gp_Pnt P = Discr.Value( i );
 330             aPnts.Append(gp_Pnt2d(P.X(), P.Y()));
 331           }
 332         }
 333         else
 334         {
 335           for( int i = Discr.NbPoints(); i > 0; i-- )
 336           {
 337             gp_Pnt P = Discr.Value( i );
 338             aPnts.Append(gp_Pnt2d(P.X(), P.Y()));
 339           }
 340         }
 341       }
 342     }
 343     NCollection_Sequence<gp_Pnt2d> aNPnts;
 344     aNPnts.Append(aPnts.First());
 345     for (int j = 1; j <= aPnts.Size() - 1; j++)
 346     {
 347       if (!aPnts(j).IsEqual(aPnts(j + 1), Precision::Confusion()))
 348         aNPnts.Append(aPnts(j + 1));
 349     }
 350
 351     //assume that the orientation of external wire & internal wires is correct
 352     //so just write all points "as-is"
 353     //External wire will be written in clockwise direction
 354     //any other wires (holes) - in anticlockwise direction
 355     for (int j = aNPnts.Size(); j >= 1; j--)
 356     {
 357       gp_XY P = gp_XY(aNPnts(j).X(), aNPnts(j).Y());
 358       theDocument->Transform(P, false);
 359       x.push_back( P.X());
 360       y.push_back( P.Y());
 361     }
 362     //first point is same as the last one => closed polygon
 363     gp_XY P = gp_XY(aNPnts.Last().X(), aNPnts.Last().Y());
 364     theDocument->Transform(P, false);
 365     x.push_back( P.X());
 366     y.push_back( P.Y());
 367
 368   }
 369
 370   aSHPObj = SHPCreateObject( SHPT_POLYGON, -1, NbWires, &anPartStart[0], NULL, x.size(), &x[0], &y[0], NULL, NULL );
 371   SHPWriteObject( theShpHandle, -1, aSHPObj );
 372   SHPDestroyObject( aSHPObj );
 373
 374 }
 375
 376 bool HYDROData_ShapeFile::Parse(SHPHandle theHandle, ShapeType theType, int& theShapeTypeOfFile)
 377 {
 378   int aShapeType;
 379   mySHPObjects.clear();
 380   SHPGetInfo( theHandle, NULL, &aShapeType, NULL, NULL );
 381   theShapeTypeOfFile = aShapeType;
 382   bool ToRead = (theType == ShapeType_Polyline &&
 383     (aShapeType == 3 || aShapeType == 13 || aShapeType == 23 || aShapeType == 5)) ||
 384    (theType == ShapeType_Polygon && aShapeType == 5);
 385   if (ToRead)
 386   {
 387     for (int i = 0; i < theHandle->nRecords; i++)
 388       mySHPObjects.push_back(SHPReadObject(theHandle, i));
 389     return true;
 390   }
 391   else
 392     return false;
 393 }
 394
 395 void HYDROData_ShapeFile::ReadSHPPolygon(Handle(HYDROData_Document) theDocument, SHPObject* anObj, int i, TopoDS_Face& F)
 396 {
 397   if (!anObj)
 398     return;
 399   TopoDS_Edge E;
 400   int nParts = anObj->nParts;
 401   Handle(Geom_Plane) aPlaneSur = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
 402
 403   BRep_Builder BB;
 404   BB.MakeFace(F);
 405
 406   NCollection_Sequence<TopoDS_Wire> allWires;
 407
 408   TopTools_SequenceOfShape aWires;
 409   for ( int i = 0 ; i < nParts ; i++ )
 410   {
 411     BRepBuilderAPI_MakeWire aBuilder;
 412     int StartIndex = anObj->panPartStart[i];
 413     int EndIndex;
 414     if (i != nParts - 1)
 415       EndIndex = anObj->panPartStart[i + 1];
 416     else
 417       EndIndex = anObj->nVertices;
 418
 419     TopoDS_Wire W;
 420     BB.MakeWire(W);
 421
 422     //First point is same as the last point
 423     int NbPnts = EndIndex - StartIndex - 1;
 424     NCollection_Array1<TopoDS_Vertex> VPoints(0, NbPnts);
 425     int j = NbPnts;
 426     for ( int k = StartIndex; k < EndIndex; k++ )
 427     {
 428       gp_Pnt P (anObj->padfX[k], anObj->padfY[k], 0);
 429       theDocument->Transform(P, true);
 430       VPoints.ChangeValue(j) = BRepLib_MakeVertex(P).Vertex();
 431       j--;
 432     }
 433
 434     for ( int k = 0; k < VPoints.Size() - 1; k++ )
 435     {
 436       gp_Pnt P1 = BRep_Tool::Pnt(VPoints(k));
 437       gp_Pnt P2 = BRep_Tool::Pnt(VPoints(k + 1));
 438       if (P1.Distance(P2) < Precision::Confusion())
 439         continue;
 440       Handle(Geom_TrimmedCurve) aTC = GC_MakeSegment(P1, P2).Value();
 441       TopoDS_Edge E;
 442       if ( k != VPoints.Size() - 2)
 443         E = BRepLib_MakeEdge(aTC, VPoints(k), VPoints(k + 1)).Edge();
 444       else
 445         //the last edge => use first and last vertices
 446         E = BRepLib_MakeEdge(aTC, VPoints.First(), VPoints.Value(VPoints.Upper() - 1)).Edge();
 447       //Add edge to wire
 448       //If SHP file is correct then the outer wire and the holes will have the correct orientations
 449       W.Closed (Standard_True);
 450       W.Orientation(TopAbs_FORWARD);
 451       BB.Add(W, E);
 452     }
 453     allWires.Append(W);
 454   }
 455
 456   int OutWIndex = -1;
 457   if (allWires.Size() > 1)
 458   {
 459     NCollection_Sequence<Bnd_Box> BBs;
 460     //try to find the largest bbox
 461     for (int i = 1; i <= allWires.Size(); i++)
 462     {
 463       TopoDS_Wire W = allWires(i);
 464       Bnd_Box BB;
 465       BRepBndLib::AddClose(W, BB);
 466       BBs.Append(BB);
 467     }
 468     for (int i = 1; i <= BBs.Size(); i++)
 469     {
 470       bool IsIn = false;
 471       for (int j = 1; j <= BBs.Size(); j++)
 472       {
 473         if (i == j)
 474           continue;
 475         Standard_Real iXmax, iXmin, iYmax, iYmin, z0, z1;
 476         Standard_Real jXmax, jXmin, jYmax, jYmin;
 477         BBs(i).Get(iXmin, iYmin, z0, iXmax, iYmax, z1);
 478         BBs(j).Get(jXmin, jYmin, z0, jXmax, jYmax, z1);
 479         if (!(iXmin > jXmin &&
 480             iYmin > jYmin &&
 481             iXmax < jXmax &&
 482             iYmax < jYmax))
 483           IsIn = true;
 484       }
 485       if (IsIn)
 486       {
 487         OutWIndex = i;
 488         break;
 489       }
 490     }
 491   }
 492   else
 493     OutWIndex = 1; //one wire => no need to check
 494
 495   for (int i = 1; i <= allWires.Size(); i++)
 496   {
 497     TopoDS_Face DF;
 498     BB.MakeFace(DF);
 499     TopoDS_Wire W = allWires(i);
 500     BB.Add(DF, W);
 501     BB.UpdateFace(DF, aPlaneSur, TopLoc_Location(), Precision::Confusion());
 502     //
 503     BRepTopAdaptor_FClass2d FClass(DF, Precision::PConfusion());
 504     if ( i == OutWIndex && FClass.PerformInfinitePoint() == TopAbs_IN)
 505       W.Reverse();
 506     if ( i != OutWIndex && FClass.PerformInfinitePoint() == TopAbs_OUT)
 507       W.Reverse();
 508     //
 509     BB.Add(F, W);
 510   }
 511
 512   //Add surface to the face
 513   BB.UpdateFace(F, aPlaneSur, TopLoc_Location(), Precision::Confusion());
 514   F.Closed(Standard_True);
 515 }
 516
 517 int HYDROData_ShapeFile::ImportPolygons(Handle(HYDROData_Document) theDocument, const QString theFileName, QStringList& thePolygonsList, TopTools_SequenceOfShape& theFaces, int& theShapeTypeOfFile)
 518 {
 519   Free();
 520   int Stat = TryOpenShapeFile(theFileName);
 521   if (Stat != 0)
 522     return Stat;
 523   myHSHP = SHPOpen( theFileName.toLatin1().data(), "rb" );
 524   if (!Parse(myHSHP, HYDROData_ShapeFile::ShapeType_Polygon, theShapeTypeOfFile))
 525     return 0;
 526   for (size_t i = 0; i < mySHPObjects.size(); i++)
 527     thePolygonsList.append("polygon_" + QString::number(i + 1));
 528
 529   TopoDS_Face aF;
 530   if (myHSHP->nShapeType == 5)
 531   {
 532 #ifdef OSD_TIMER
 533     OSD_Timer timer;
 534     timer.Start();
 535 #endif
 536     for (size_t i = 0; i < mySHPObjects.size(); i++)
 537     {
 538        ReadSHPPolygon(theDocument, mySHPObjects[i], i, aF);
 539        theFaces.Append(aF);
 540     }
 541 #ifdef OSD_TIMER
 542     timer.Stop();
 543     timer.Show();
 544 #endif
 545     return 1;
 546   }
 547   else
 548     return 0;
 549 }
 550
 551 void HYDROData_ShapeFile::Free()
 552 {
 553   for (size_t i = 0; i < mySHPObjects.size(); i++ )
 554     free (mySHPObjects[i]);
 555
 556   mySHPObjects.clear();
 557   if (myHSHP != NULL)
 558   {
 559     SHPClose(myHSHP);
 560     myHSHP = NULL;
 561   }
 562 }
 563
 564
 565 void HYDROData_ShapeFile::ReadSHPPolyXY(Handle(HYDROData_Document) theDocument, SHPObject* anObj, QString theFileName,
 566   int theInd, NCollection_Sequence<Handle(HYDROData_Entity)>& theEntities, bool bReadAsPolyline)
 567 {
 568
 569   Handle(HYDROData_PolylineXY) aPolylineXY = Handle(HYDROData_PolylineXY)::DownCast( theDocument->CreateObject( KIND_POLYLINEXY ) );
 570
 571   int nParts = anObj->nParts;
 572   for ( int i = 0 ; i < nParts ; i++ )
 573   {
 574     int StartIndex = anObj->panPartStart[i];
 575     int EndIndex;
 576     if (i != nParts - 1)
 577       EndIndex = anObj->panPartStart[i + 1];
 578     else
 579       EndIndex = anObj->nVertices;
 580
 581     HYDROData_PolylineXY::SectionType aSectType = HYDROData_PolylineXY::SECTION_POLYLINE;
 582     if (!bReadAsPolyline)
 583     {
 584       bool IsClosed = false;
 585       if (anObj->padfX[StartIndex] == anObj->padfX[EndIndex - 1] &&
 586         anObj->padfY[StartIndex] == anObj->padfY[EndIndex - 1] )
 587       {
 588         IsClosed = true;
 589         aPolylineXY->AddSection( TCollection_AsciiString( ("poly_section_" + QString::number(i)).data()->toLatin1()), aSectType, true);
 590       }
 591       else
 592         aPolylineXY->AddSection( TCollection_AsciiString( ("poly_section_" + QString::number(i)).data()->toLatin1()), aSectType, false);
 593
 594       if (IsClosed)
 595         EndIndex--;
 596     }
 597     else
 598     {
 599       //polygon; contours always closed
 600       aPolylineXY->AddSection( TCollection_AsciiString( ("poly_section_" + QString::number(i)).data()->toLatin1()), aSectType, true);
 601       EndIndex--;
 602     }
 603     for ( int k = StartIndex; k < EndIndex ; k++ )
 604     {
 605       HYDROData_PolylineXY::Point aSectPoint = gp_XY(anObj->padfX[k], anObj->padfY[k]);
 606       theDocument->Transform(aSectPoint, true);
 607       aPolylineXY->AddPoint( i, aSectPoint );
 608     }
 609
 610   }
 611
 612   aPolylineXY->SetWireColor( HYDROData_PolylineXY::DefaultWireColor() );
 613   aPolylineXY->SetName( theFileName + "_PolyXY_" + QString::number(theInd) );
 614
 615   aPolylineXY->Update();
 616   theEntities.Append(aPolylineXY);
 617
 618 }
 619
 620 void HYDROData_ShapeFile::ReadSHPPoly3D(Handle(HYDROData_Document) theDocument, SHPObject* anObj, QString theFileName,
 621   int theInd, NCollection_Sequence<Handle(HYDROData_Entity)>& theEntities)
 622 {
 623   Handle(HYDROData_PolylineXY) aPolylineXY = Handle(HYDROData_PolylineXY)::DownCast( theDocument->CreateObject( KIND_POLYLINEXY ) );
 624
 625   Handle(HYDROData_Polyline3D) aPolylineObj = Handle(HYDROData_Polyline3D)::DownCast( theDocument->CreateObject( KIND_POLYLINE ) );
 626
 627   Handle(HYDROData_Bathymetry) aBath = Handle(HYDROData_Bathymetry)::DownCast( theDocument->CreateObject( KIND_BATHYMETRY ) );
 628   HYDROData_Bathymetry::AltitudePoints aAPoints;
 629
 630   int nParts = anObj->nParts;
 631   for ( int i = 0 ; i < nParts ; i++ )
 632   {
 633     //bool aSectClosure = true;
 634     int StartIndex = anObj->panPartStart[i];
 635     int EndIndex;
 636     if (i != nParts - 1)
 637       EndIndex = anObj->panPartStart[i + 1];
 638     else
 639       EndIndex = anObj->nVertices;
 640
 641     bool IsClosed = false;
 642     HYDROData_PolylineXY::SectionType aSectType = HYDROData_PolylineXY::SECTION_POLYLINE;
 643     if (anObj->padfX[StartIndex] == anObj->padfX[EndIndex - 1] &&
 644         anObj->padfY[StartIndex] == anObj->padfY[EndIndex - 1] &&
 645         anObj->padfZ[StartIndex] == anObj->padfZ[EndIndex - 1])
 646     {
 647       IsClosed = true;
 648       aPolylineXY->AddSection( TCollection_AsciiString( ("poly_section_" + QString::number(i)).data()->toLatin1()), aSectType, true );
 649     }
 650     else
 651       aPolylineXY->AddSection( TCollection_AsciiString( ("poly_section_" + QString::number(i)).data()->toLatin1()), aSectType, false  );
 652
 653     if (IsClosed)
 654       EndIndex--;
 655     for ( int k = StartIndex ; k < EndIndex ; k++ )
 656     {
 657       HYDROData_PolylineXY::Point aSectPoint = gp_XY(anObj->padfX[k], anObj->padfY[k]);
 658       theDocument->Transform(aSectPoint, true);
 659       aPolylineXY->AddPoint( i, aSectPoint );
 660       HYDROData_Bathymetry::AltitudePoint p;
 661       p.X = aSectPoint.X();
 662       p.Y = aSectPoint.Y();
 663       p.Z = anObj->padfZ[k];
 664       aAPoints.push_back(p);
 665     }
 666   }
 667
 668
 669   QString aBathName = theFileName + "_bath_" + QString::number(theInd);
 670   QString aPolyXYName = theFileName + "_polyXY_" + QString::number(theInd);
 671   QString aPoly3DName = theFileName + "_poly3D_" + QString::number(theInd);
 672
 673   aPolylineXY->SetName( aPolyXYName );
 674   aPolylineXY->SetWireColor(HYDROData_PolylineXY::DefaultWireColor());
 675   aPolylineXY->Update();
 676
 677   aBath->SetAltitudePoints(aAPoints);
 678   aBath->SetName( aBathName );
 679
 680   aPolylineObj->SetPolylineXY (aPolylineXY, false);
 681   aPolylineObj->SetAltitudeObject(aBath);
 682
 683   aPolylineObj->SetBorderColor( aPolylineObj->DefaultBorderColor() );
 684   aPolylineObj->SetName( aPoly3DName );
 685
 686   aPolylineObj->Update();
 687   theEntities.Append(aPolylineXY);
 688   theEntities.Append(aPolylineObj);
 689
 690 }
 691
 692 void HYDROData_ShapeFile::GetFreeIndices(std::vector<int>& theAllowedIndexes, QString strName, size_t theObjsSize,
 693                                          QStringList theExistingNames, QString theBaseFileName)
 694 {
 695   int anInd = 0;
 696   for (;theAllowedIndexes.size() < theObjsSize;)
 697   {
 698     if (!theExistingNames.contains(theBaseFileName + strName + QString::number(anInd)))
 699     {
 700       theAllowedIndexes.push_back(anInd);
 701       anInd++;
 702     }
 703     else
 704       anInd++;
 705   }
 706 }
 707
 708 int HYDROData_ShapeFile::ImportPolylines(Handle(HYDROData_Document) theDocument, const QString& theFileName,
 709   NCollection_Sequence<Handle(HYDROData_Entity)>& theEntities, int& theShapeTypeOfFile)
 710 {
 711   //Free();
 712   int aStat = TryOpenShapeFile(theFileName);
 713   if (aStat != 0)
 714     return aStat;
 715
 716   HYDROData_Iterator anIter( theDocument );
 717   int anInd = 0;
 718   QStringList anExistingNames;
 719   std::vector<int> anAllowedIndexes;
 720   for( ; anIter.More(); anIter.Next() )
 721     anExistingNames.push_back(anIter.Current()->GetName());
 722
 723   SHPHandle aHSHP;
 724   aHSHP = SHPOpen( theFileName.toLatin1().data(), "rb" );
 725
 726   QFileInfo aFileInfo(theFileName);
 727   QString aBaseFileName = aFileInfo.baseName();
 728
 729   if (!Parse(aHSHP, HYDROData_ShapeFile::ShapeType_Polyline, theShapeTypeOfFile))
 730     return 0;
 731   if (aHSHP->nShapeType == 3 || aHSHP->nShapeType == 23)
 732   {
 733     size_t anObjsSize = mySHPObjects.size();
 734     GetFreeIndices(anAllowedIndexes, "_PolyXY_", anObjsSize, anExistingNames, aBaseFileName);
 735
 736     for (size_t i = 0; i < mySHPObjects.size(); i++ )
 737     {
 738       ReadSHPPolyXY(theDocument, mySHPObjects[i], aBaseFileName, anAllowedIndexes[i], theEntities, false);
 739     }
 740     aStat = 1;
 741   }
 742   else if (aHSHP->nShapeType == 13)
 743   {
 744     anInd = 0;
 745     for (;anAllowedIndexes.size() < mySHPObjects.size();)
 746     {
 747       if (!anExistingNames.contains(aBaseFileName + "_PolyXY_" + QString::number(anInd)) &&
 748           !anExistingNames.contains(aBaseFileName + "_Poly3D_" + QString::number(anInd)) &&
 749           !anExistingNames.contains(aBaseFileName + "_Bath_" + QString::number(anInd)))
 750       {
 751         anAllowedIndexes.push_back(anInd);
 752         anInd++;
 753       }
 754       else
 755         anInd++;
 756     }
 757     for (size_t i = 0; i < mySHPObjects.size(); i++ )
 758     {
 759       ReadSHPPoly3D(theDocument, mySHPObjects[i], aBaseFileName, anAllowedIndexes[i], theEntities);
 760     }
 761     aStat = 1;
 762   }
 763   else if (aHSHP->nShapeType == 5)
 764   {
 765     //import polygon's contours as polylines
 766     size_t anObjsSize = mySHPObjects.size();
 767     GetFreeIndices(anAllowedIndexes, "_PolyXY_", anObjsSize, anExistingNames, aBaseFileName);
 768
 769     for (size_t i = 0; i < mySHPObjects.size(); i++ )
 770     {
 771       ReadSHPPolyXY(theDocument, mySHPObjects[i], aBaseFileName, anAllowedIndexes[i], theEntities, true);
 772     }
 773     aStat = 2;
 774   }
 775   else
 776   {
 777     aStat = 0;
 778   }
 779
 780   for (size_t i = 0; i < mySHPObjects.size(); i++ )
 781     free (mySHPObjects[i]);
 782
 783   mySHPObjects.clear();
 784   SHPClose(aHSHP);
 785   return aStat;
 786 }
 787
 788 QString HYDROData_ShapeFile::GetShapeTypeName(int theType)
 789 {
 790   switch (theType)
 791   {
 792     case 0:
 793       return "null shape";
 794     case 1:
 795       return "point (unsupported by HYDRO)";
 796     case 3:
 797       return "arc/polyline (supported by HYDRO)";
 798     case 5:
 799       return "polygon (supported by HYDRO)";
 800     case 8:
 801       return "multipoint (unsupported by HYDRO)";
 802     case 11:
 803       return "pointZ (unsupported by HYDRO)";
 804     case 13:
 805       return "arcZ/polyline (supported by HYDRO)";
 806     case 15:
 807       return "polygonZ (unsupported by HYDRO)";
 808     case 18:
 809       return "multipointZ (unsupported by HYDRO)";
 810     case 21:
 811       return "pointM (unsupported by HYDRO)";
 812     case 23:
 813       return "arcM/polyline (supported by HYDRO)";
 814     case 25:
 815       return "polygonM (unsupported by HYDRO)";
 816     case 28:
 817       return "multipointM (unsupported by HYDRO)";
 818     case 31:
 819       return "multipatch (unsupported by HYDRO)";
 820     default:
 821       return "unknown";
 822   }
 823 }
 824
 825 int HYDROData_ShapeFile::TryOpenShapeFile(QString theFileName)
 826 {
 827   QString aSHPfile = theFileName.simplified();
 828   QString aSHXfile = theFileName.simplified().replace( theFileName.simplified().size() - 4, 4, ".shx");
 829
 830   QString anExt = theFileName.split('.', QString::SkipEmptyParts).back();
 831   if (anExt.toLower() != "shp")
 832     return -3;
 833
 834   FILE* pFileSHP = NULL;
 835   pFileSHP = fopen (aSHPfile.toLatin1().data(), "r");
 836   FILE* pFileSHX = NULL;
 837   pFileSHX = fopen (aSHXfile.toLatin1().data(), "r");
 838
 839   if (pFileSHP == NULL || pFileSHX == NULL)
 840   {
 841     if (pFileSHP == NULL)
 842       return -1;
 843     if (pFileSHX == NULL)
 844       return -2;
 845   }
 846
 847   fclose (pFileSHP);
 848   fclose (pFileSHX);
 849   return 0;
 850 }
 851
 852
 853 bool HYDROData_ShapeFile::CheckDBFFileExisting(const QString& theSHPFilePath, QString& thePathToDBFFile)
 854 {
 855   QString aSHPfile = theSHPFilePath.simplified();
 856   QString aDBFfile = theSHPFilePath.simplified().replace( theSHPFilePath.simplified().size() - 4, 4, ".dbf");
 857   FILE* pFileDBF = NULL;
 858   pFileDBF = fopen (aDBFfile.toLatin1().data(), "r");
 859
 860   if (pFileDBF == NULL)
 861   {
 862     return false;
 863   }
 864
 865   fclose (pFileDBF);
 866   thePathToDBFFile = aDBFfile;
 867   return true;
 868 }
 869
 870
 871 bool HYDROData_ShapeFile::DBF_OpenDBF(const QString& thePathToDBFFile)
 872 {
 873   myHDBF = DBFOpen( thePathToDBFFile.toLatin1().data(), "r" );
 874   if(myHDBF != NULL)
 875     return true;
 876   else
 877     return false;
 878 }
 879
 880 int HYDROData_ShapeFile::DBF_GetNbFields()
 881 {
 882   if(myHDBF == NULL)
 883     return 0;
 884   return DBFGetFieldCount(myHDBF);
 885 }
 886
 887 void HYDROData_ShapeFile::DBF_CloseDBF()
 888 {
 889   if(myHDBF != NULL)
 890      DBFClose( myHDBF );
 891 }
 892
 893 QStringList HYDROData_ShapeFile::DBF_GetFieldList()
 894 {
 895   QStringList FieldList;
 896   int   nWidth, nDecimals;
 897   char chField[12];
 898
 899   for( int i = 0; i < DBFGetFieldCount(myHDBF); i++ )
 900   {
 901       DBFFieldType eType;
 902       eType = DBFGetFieldInfo( myHDBF, i, chField, &nWidth, &nDecimals );
 903       FieldList.append(QString(chField));
 904   }
 905
 906   return FieldList;
 907 }
 908
 909 void HYDROData_ShapeFile::DBF_GetFieldTypeList(std::vector<DBF_FieldType>& FTVect)
 910 {
 911   int   nWidth, nDecimals;
 912   char chField[12];
 913   DBF_FieldType FT;
 914   for( int i = 0; i < DBFGetFieldCount(myHDBF); i++ )
 915   {
 916     DBFFieldType eType;
 917     eType = DBFGetFieldInfo( myHDBF, i, chField, &nWidth, &nDecimals );
 918     if( eType == FTString )
 919       FT = DBF_FieldType_String;
 920     else if( eType == FTInteger )
 921       FT = DBF_FieldType_Integer;
 922     else if( eType == FTDouble )
 923       FT = DBF_FieldType_Double;
 924     else if( eType == FTInvalid )
 925       FT = DBF_FieldType_Invalid;
 926
 927     FTVect.push_back(FT);
 928   }
 929
 930 }
 931
 932 int HYDROData_ShapeFile::DBF_GetNbRecords()
 933 {
 934   if(myHDBF == NULL)
 935     return 0;
 936   return DBFGetRecordCount(myHDBF);
 937 }
 938
 939 void HYDROData_ShapeFile::DBF_GetAttributeList(int theIndexOfField, std::vector<DBF_AttrValue>& theAttrV)
 940 {
 941   int nWidth, nDecimals;
 942   char chField[12];
 943
 944   for( int i = 0; i < DBFGetRecordCount(myHDBF); i++ )
 945   {
 946     DBFFieldType eType;
 947     DBF_AttrValue anAttr;
 948     eType = DBFGetFieldInfo( myHDBF, theIndexOfField, chField, &nWidth, &nDecimals );
 949
 950     if( DBFIsAttributeNULL( myHDBF, i, theIndexOfField ) )
 951     {
 952       anAttr.myIsNull = true;
 953       DBF_FieldType FT = DBF_FieldType_None;
 954       if( eType == FTString )
 955         FT = DBF_FieldType_String;
 956       else if( eType == FTInteger )
 957         FT = DBF_FieldType_Integer;
 958       else if( eType == FTDouble )
 959         FT = DBF_FieldType_Double;
 960       else if( eType == FTInvalid )
 961         FT = DBF_FieldType_Invalid;
 962       anAttr.myFieldType = FT;
 963     }
 964     else
 965     {
 966       switch( eType )
 967       {
 968         case FTString:
 969         {
 970           const char* chAttr = DBFReadStringAttribute( myHDBF, i, theIndexOfField );
 971           anAttr.myIsNull = false;
 972           anAttr.myFieldType = DBF_FieldType_String;
 973           anAttr.myRawValue = chAttr;
 974           anAttr.myStrVal = QString(chAttr);
 975           break;
 976         }
 977
 978         case FTInteger:
 979         {
 980           int iAttr = DBFReadIntegerAttribute( myHDBF, i, theIndexOfField );
 981           anAttr.myIsNull = false;
 982           anAttr.myFieldType = DBF_FieldType_Integer;
 983           anAttr.myRawValue = DBFReadStringAttribute( myHDBF, i, theIndexOfField );
 984           anAttr.myIntVal = iAttr;
 985           break;
 986         }
 987
 988         case FTDouble:
 989         {
 990           double dAttr = DBFReadDoubleAttribute( myHDBF, i, theIndexOfField );
 991           anAttr.myIsNull = false;
 992           anAttr.myFieldType = DBF_FieldType_Double;
 993           anAttr.myRawValue = DBFReadStringAttribute( myHDBF, i, theIndexOfField );
 994           anAttr.myDoubleVal = dAttr;
 995           break;
 996         }
 997         default:
 998           break;
 999       }
1000     }
1001     theAttrV.push_back(anAttr);
1002   }
1003
1004 }
1005
1006 bool HYDROData_ShapeFile::DBF_WriteFieldAndValues(const QString& theFileName, const QString& theFieldName, DBF_FieldType theType, const std::vector<DBF_AttrValue>& theAttrV, bool bUseRawValue)
1007 {
1008   // Check that given field type is equal to field types of attributes values
1009   for (size_t i = 0; i < theAttrV.size(); i++)
1010   {
1011     if (theAttrV[i].myFieldType != theType)
1012       return false;
1013   }
1014
1015   DBFHandle     hDBF;
1016   hDBF = DBFCreate( theFileName.toStdString().c_str() );
1017   if( hDBF == NULL )
1018           return false;
1019
1020   if (theType != DBF_FieldType_String && theType != DBF_FieldType_Integer && theType != DBF_FieldType_Double)
1021   {
1022     DBFClose( hDBF );
1023     return false; //cant handle any other cases
1024   }
1025
1026   int nWidth = 20;
1027   switch( theType )
1028   {
1029     case DBF_FieldType_String:
1030     {
1031       DBFAddField (hDBF, theFieldName.toStdString().c_str(), FTString, nWidth, 0);
1032       break;
1033     }
1034
1035     case DBF_FieldType_Integer:
1036     {
1037       DBFAddField (hDBF, theFieldName.toStdString().c_str(), FTInteger, nWidth, 0);
1038       break;
1039     }
1040
1041     case DBF_FieldType_Double:
1042     {
1043       DBFAddField (hDBF, theFieldName.toStdString().c_str(), FTDouble, nWidth, 0);
1044       break;
1045     }
1046     default:
1047       break;
1048   }
1049
1050   if (DBFGetFieldCount( hDBF ) != 1)
1051   {
1052     DBFClose( hDBF );
1053     return false;
1054   }
1055   if (DBFGetRecordCount( hDBF ) != 0)
1056   {
1057     DBFClose( hDBF );
1058     return false;
1059   }
1060   int stat = -1;
1061
1062   if (bUseRawValue)
1063   {
1064     for (size_t i = 0; i < theAttrV.size(); i++)
1065     {
1066       if (!theAttrV[i].myIsNull)
1067         stat = DBFWriteStringAttribute(hDBF, (int)i, 0, theAttrV[i].myRawValue.c_str());
1068       else
1069         stat = DBFWriteNULLAttribute(hDBF, (int)i, 0 );
1070
1071       if (stat != 1)
1072       {
1073         DBFClose( hDBF );
1074         return false;
1075       }
1076     }
1077   }
1078   else
1079   {
1080     for (size_t i = 0; i < theAttrV.size(); i++)
1081     {
1082       if (!theAttrV[i].myIsNull)
1083         switch( theType )
1084         {
1085           case DBF_FieldType_String:
1086           {
1087             stat = DBFWriteStringAttribute(hDBF, (int)i, 0, theAttrV[i].myStrVal.toStdString().c_str());
1088             break;
1089           }
1090
1091           case DBF_FieldType_Integer:
1092           {
1093             stat = DBFWriteIntegerAttribute(hDBF, (int)i, 0, theAttrV[i].myIntVal);
1094             break;
1095           }
1096
1097           case DBF_FieldType_Double:
1098           {
1099             stat = DBFWriteDoubleAttribute(hDBF, (int)i, 0, theAttrV[i].myDoubleVal);
1100             break;
1101           }
1102           default:
1103             break;
1104         }
1105       else
1106         stat = DBFWriteNULLAttribute(hDBF, (int)i, 0 );
1107
1108       if (stat != 1)
1109       {
1110         DBFClose( hDBF );
1111         return false;
1112       }
1113     }
1114   }
1115
1116   DBFClose( hDBF );
1117   return true;
1118
1119 }
1120