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, ImportShapeType theShapeTypesToImport)
 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
 732   int sh_type = aHSHP->nShapeType;
 733
 734   if ((theShapeTypesToImport == HYDROData_ShapeFile::ImportShapeType_All ||
 735     theShapeTypesToImport == HYDROData_ShapeFile::ImportShapeType_Polyline)
 736     && (sh_type == 3 || sh_type == 23))
 737   {
 738     size_t anObjsSize = mySHPObjects.size();
 739     GetFreeIndices(anAllowedIndexes, "_PolyXY_", anObjsSize, anExistingNames, aBaseFileName);
 740
 741     for (size_t i = 0; i < mySHPObjects.size(); i++ )
 742     {
 743       ReadSHPPolyXY(theDocument, mySHPObjects[i], aBaseFileName, anAllowedIndexes[i], theEntities, false);
 744     }
 745     aStat = 1;
 746   }
 747   else if ((theShapeTypesToImport == HYDROData_ShapeFile::ImportShapeType_All ||
 748     theShapeTypesToImport == HYDROData_ShapeFile::ImportShapeType_Polyline3D) && sh_type == 13)
 749   {
 750     anInd = 0;
 751     for (;anAllowedIndexes.size() < mySHPObjects.size();)
 752     {
 753       if (!anExistingNames.contains(aBaseFileName + "_PolyXY_" + QString::number(anInd)) &&
 754           !anExistingNames.contains(aBaseFileName + "_Poly3D_" + QString::number(anInd)) &&
 755           !anExistingNames.contains(aBaseFileName + "_Bath_" + QString::number(anInd)))
 756       {
 757         anAllowedIndexes.push_back(anInd);
 758         anInd++;
 759       }
 760       else
 761         anInd++;
 762     }
 763     for (size_t i = 0; i < mySHPObjects.size(); i++ )
 764     {
 765       ReadSHPPoly3D(theDocument, mySHPObjects[i], aBaseFileName, anAllowedIndexes[i], theEntities);
 766     }
 767     aStat = 1;
 768   }
 769   else if ((theShapeTypesToImport == HYDROData_ShapeFile::ImportShapeType_All ||
 770     theShapeTypesToImport == HYDROData_ShapeFile::ImportShapeType_Polygon) && sh_type == 5)
 771   {
 772     //import polygon's contours as polylines
 773     size_t anObjsSize = mySHPObjects.size();
 774     GetFreeIndices(anAllowedIndexes, "_PolyXY_", anObjsSize, anExistingNames, aBaseFileName);
 775
 776     for (size_t i = 0; i < mySHPObjects.size(); i++ )
 777     {
 778       ReadSHPPolyXY(theDocument, mySHPObjects[i], aBaseFileName, anAllowedIndexes[i], theEntities, true);
 779     }
 780     aStat = 2;
 781   }
 782   else
 783   {
 784     aStat = 0;
 785   }
 786
 787   for (size_t i = 0; i < mySHPObjects.size(); i++ )
 788     free (mySHPObjects[i]);
 789
 790   mySHPObjects.clear();
 791   SHPClose(aHSHP);
 792   return aStat;
 793 }
 794
 795 QString HYDROData_ShapeFile::GetShapeTypeName(int theType)
 796 {
 797   switch (theType)
 798   {
 799     case 0:
 800       return "null shape";
 801     case 1:
 802       return "point (unsupported by HYDRO)";
 803     case 3:
 804       return "arc/polyline (supported by HYDRO)";
 805     case 5:
 806       return "polygon (supported by HYDRO)";
 807     case 8:
 808       return "multipoint (unsupported by HYDRO)";
 809     case 11:
 810       return "pointZ (unsupported by HYDRO)";
 811     case 13:
 812       return "arcZ/polyline (supported by HYDRO)";
 813     case 15:
 814       return "polygonZ (unsupported by HYDRO)";
 815     case 18:
 816       return "multipointZ (unsupported by HYDRO)";
 817     case 21:
 818       return "pointM (unsupported by HYDRO)";
 819     case 23:
 820       return "arcM/polyline (supported by HYDRO)";
 821     case 25:
 822       return "polygonM (unsupported by HYDRO)";
 823     case 28:
 824       return "multipointM (unsupported by HYDRO)";
 825     case 31:
 826       return "multipatch (unsupported by HYDRO)";
 827     default:
 828       return "unknown";
 829   }
 830 }
 831
 832 int HYDROData_ShapeFile::TryOpenShapeFile(QString theFileName)
 833 {
 834   QString aSHPfile = theFileName.simplified();
 835   QString aSHXfile = theFileName.simplified().replace( theFileName.simplified().size() - 4, 4, ".shx");
 836
 837   QString anExt = theFileName.split('.', QString::SkipEmptyParts).back();
 838   if (anExt.toLower() != "shp")
 839     return -3;
 840
 841   FILE* pFileSHP = NULL;
 842   pFileSHP = fopen (aSHPfile.toLatin1().data(), "r");
 843   FILE* pFileSHX = NULL;
 844   pFileSHX = fopen (aSHXfile.toLatin1().data(), "r");
 845
 846   if (pFileSHP == NULL || pFileSHX == NULL)
 847   {
 848     if (pFileSHP == NULL)
 849       return -1;
 850     if (pFileSHX == NULL)
 851       return -2;
 852   }
 853
 854   fclose (pFileSHP);
 855   fclose (pFileSHX);
 856   return 0;
 857 }
 858
 859
 860 bool HYDROData_ShapeFile::CheckDBFFileExisting(const QString& theSHPFilePath, QString& thePathToDBFFile)
 861 {
 862   QString aSHPfile = theSHPFilePath.simplified();
 863   QString aDBFfile = theSHPFilePath.simplified().replace( theSHPFilePath.simplified().size() - 4, 4, ".dbf");
 864   FILE* pFileDBF = NULL;
 865   pFileDBF = fopen (aDBFfile.toLatin1().data(), "r");
 866
 867   if (pFileDBF == NULL)
 868   {
 869     return false;
 870   }
 871
 872   fclose (pFileDBF);
 873   thePathToDBFFile = aDBFfile;
 874   return true;
 875 }
 876
 877
 878 bool HYDROData_ShapeFile::DBF_OpenDBF(const QString& thePathToDBFFile)
 879 {
 880   myHDBF = DBFOpen( thePathToDBFFile.toLatin1().data(), "r" );
 881   if(myHDBF != NULL)
 882     return true;
 883   else
 884     return false;
 885 }
 886
 887 int HYDROData_ShapeFile::DBF_GetNbFields()
 888 {
 889   if(myHDBF == NULL)
 890     return 0;
 891   return DBFGetFieldCount(myHDBF);
 892 }
 893
 894 void HYDROData_ShapeFile::DBF_CloseDBF()
 895 {
 896   if(myHDBF != NULL)
 897      DBFClose( myHDBF );
 898 }
 899
 900 QStringList HYDROData_ShapeFile::DBF_GetFieldList()
 901 {
 902   QStringList FieldList;
 903   int   nWidth, nDecimals;
 904   char chField[12];
 905
 906   for( int i = 0; i < DBFGetFieldCount(myHDBF); i++ )
 907   {
 908       DBFFieldType eType;
 909       eType = DBFGetFieldInfo( myHDBF, i, chField, &nWidth, &nDecimals );
 910       FieldList.append(QString(chField));
 911   }
 912
 913   return FieldList;
 914 }
 915
 916 void HYDROData_ShapeFile::DBF_GetFieldTypeList(std::vector<DBF_FieldType>& FTVect)
 917 {
 918   int   nWidth, nDecimals;
 919   char chField[12];
 920   DBF_FieldType FT;
 921   for( int i = 0; i < DBFGetFieldCount(myHDBF); i++ )
 922   {
 923     DBFFieldType eType;
 924     eType = DBFGetFieldInfo( myHDBF, i, chField, &nWidth, &nDecimals );
 925     if( eType == FTString )
 926       FT = DBF_FieldType_String;
 927     else if( eType == FTInteger )
 928       FT = DBF_FieldType_Integer;
 929     else if( eType == FTDouble )
 930       FT = DBF_FieldType_Double;
 931     else if( eType == FTInvalid )
 932       FT = DBF_FieldType_Invalid;
 933
 934     FTVect.push_back(FT);
 935   }
 936
 937 }
 938
 939 int HYDROData_ShapeFile::DBF_GetNbRecords()
 940 {
 941   if(myHDBF == NULL)
 942     return 0;
 943   return DBFGetRecordCount(myHDBF);
 944 }
 945
 946 void HYDROData_ShapeFile::DBF_GetAttributeList(int theIndexOfField, std::vector<DBF_AttrValue>& theAttrV)
 947 {
 948   int nWidth, nDecimals;
 949   char chField[12];
 950
 951   for( int i = 0; i < DBFGetRecordCount(myHDBF); i++ )
 952   {
 953     DBFFieldType eType;
 954     DBF_AttrValue anAttr;
 955     eType = DBFGetFieldInfo( myHDBF, theIndexOfField, chField, &nWidth, &nDecimals );
 956
 957     if( DBFIsAttributeNULL( myHDBF, i, theIndexOfField ) )
 958     {
 959       anAttr.myIsNull = true;
 960       DBF_FieldType FT = DBF_FieldType_None;
 961       if( eType == FTString )
 962         FT = DBF_FieldType_String;
 963       else if( eType == FTInteger )
 964         FT = DBF_FieldType_Integer;
 965       else if( eType == FTDouble )
 966         FT = DBF_FieldType_Double;
 967       else if( eType == FTInvalid )
 968         FT = DBF_FieldType_Invalid;
 969       anAttr.myFieldType = FT;
 970     }
 971     else
 972     {
 973       switch( eType )
 974       {
 975         case FTString:
 976         {
 977           const char* chAttr = DBFReadStringAttribute( myHDBF, i, theIndexOfField );
 978           anAttr.myIsNull = false;
 979           anAttr.myFieldType = DBF_FieldType_String;
 980           anAttr.myRawValue = chAttr;
 981           anAttr.myStrVal = QString(chAttr);
 982           break;
 983         }
 984
 985         case FTInteger:
 986         {
 987           int iAttr = DBFReadIntegerAttribute( myHDBF, i, theIndexOfField );
 988           anAttr.myIsNull = false;
 989           anAttr.myFieldType = DBF_FieldType_Integer;
 990           anAttr.myRawValue = DBFReadStringAttribute( myHDBF, i, theIndexOfField );
 991           anAttr.myIntVal = iAttr;
 992           break;
 993         }
 994
 995         case FTDouble:
 996         {
 997           double dAttr = DBFReadDoubleAttribute( myHDBF, i, theIndexOfField );
 998           anAttr.myIsNull = false;
 999           anAttr.myFieldType = DBF_FieldType_Double;
1000           anAttr.myRawValue = DBFReadStringAttribute( myHDBF, i, theIndexOfField );
1001           anAttr.myDoubleVal = dAttr;
1002           break;
1003         }
1004         default:
1005           break;
1006       }
1007     }
1008     theAttrV.push_back(anAttr);
1009   }
1010
1011 }
1012
1013 bool HYDROData_ShapeFile::DBF_WriteFieldAndValues(const QString& theFileName, const QString& theFieldName, DBF_FieldType theType, const std::vector<DBF_AttrValue>& theAttrV, bool bUseRawValue)
1014 {
1015   // Check that given field type is equal to field types of attributes values
1016   for (size_t i = 0; i < theAttrV.size(); i++)
1017   {
1018     if (theAttrV[i].myFieldType != theType)
1019       return false;
1020   }
1021
1022   DBFHandle     hDBF;
1023   hDBF = DBFCreate( theFileName.toStdString().c_str() );
1024   if( hDBF == NULL )
1025           return false;
1026
1027   if (theType != DBF_FieldType_String && theType != DBF_FieldType_Integer && theType != DBF_FieldType_Double)
1028   {
1029     DBFClose( hDBF );
1030     return false; //cant handle any other cases
1031   }
1032
1033   int nWidth = 20;
1034   switch( theType )
1035   {
1036     case DBF_FieldType_String:
1037     {
1038       DBFAddField (hDBF, theFieldName.toStdString().c_str(), FTString, nWidth, 0);
1039       break;
1040     }
1041
1042     case DBF_FieldType_Integer:
1043     {
1044       DBFAddField (hDBF, theFieldName.toStdString().c_str(), FTInteger, nWidth, 0);
1045       break;
1046     }
1047
1048     case DBF_FieldType_Double:
1049     {
1050       DBFAddField (hDBF, theFieldName.toStdString().c_str(), FTDouble, nWidth, 0);
1051       break;
1052     }
1053     default:
1054       break;
1055   }
1056
1057   if (DBFGetFieldCount( hDBF ) != 1)
1058   {
1059     DBFClose( hDBF );
1060     return false;
1061   }
1062   if (DBFGetRecordCount( hDBF ) != 0)
1063   {
1064     DBFClose( hDBF );
1065     return false;
1066   }
1067   int stat = -1;
1068
1069   if (bUseRawValue)
1070   {
1071     for (size_t i = 0; i < theAttrV.size(); i++)
1072     {
1073       if (!theAttrV[i].myIsNull)
1074         stat = DBFWriteStringAttribute(hDBF, (int)i, 0, theAttrV[i].myRawValue.c_str());
1075       else
1076         stat = DBFWriteNULLAttribute(hDBF, (int)i, 0 );
1077
1078       if (stat != 1)
1079       {
1080         DBFClose( hDBF );
1081         return false;
1082       }
1083     }
1084   }
1085   else
1086   {
1087     for (size_t i = 0; i < theAttrV.size(); i++)
1088     {
1089       if (!theAttrV[i].myIsNull)
1090         switch( theType )
1091         {
1092           case DBF_FieldType_String:
1093           {
1094             stat = DBFWriteStringAttribute(hDBF, (int)i, 0, theAttrV[i].myStrVal.toStdString().c_str());
1095             break;
1096           }
1097
1098           case DBF_FieldType_Integer:
1099           {
1100             stat = DBFWriteIntegerAttribute(hDBF, (int)i, 0, theAttrV[i].myIntVal);
1101             break;
1102           }
1103
1104           case DBF_FieldType_Double:
1105           {
1106             stat = DBFWriteDoubleAttribute(hDBF, (int)i, 0, theAttrV[i].myDoubleVal);
1107             break;
1108           }
1109           default:
1110             break;
1111         }
1112       else
1113         stat = DBFWriteNULLAttribute(hDBF, (int)i, 0 );
1114
1115       if (stat != 1)
1116       {
1117         DBFClose( hDBF );
1118         return false;
1119       }
1120     }
1121   }
1122
1123   DBFClose( hDBF );
1124   return true;
1125
1126 }
1127