[vtkusers] Re2: How can I add ESRI Shape Files supporting to VTK?

[Alexey Yutkin]

> Have you read the shapefile white paper on the ESRI site?
Yes.



My C++ way:

  I am download latest source code from Frank Warmerdam Shapelib library and
"ESRI Shapefile Technical Description" from ESRI web site.
 I put source code from this library into VTK/Utilities/shapelib and make
CMakeLists.txt in this directory.
-----------------------------------------------------------
 File VTK/Utilities/shapelib/CMakeLists.txt:
-----------------------------------------------------------

PROJECT(VTKSHP)

# source files for shp
SOURCE_FILES(SHP_SRCS
dbfopen.c    shpopen.c
)

CONFIGURE_FILE(
${VTKSHP_SOURCE_DIR}/.NoDartCoverage
${VTKSHP_BINARY_DIR}/.NoDartCoverage)

IF(WIN32)
  IF(BUILD_SHARED_LIBS)
    ADD_DEFINITIONS(-DSHAPELIB_DLLEXPORT)
  ELSE(BUILD_SHARED_LIBS)
    ADD_DEFINITIONS()
  ENDIF(BUILD_SHARED_LIBS)
ENDIF(WIN32)

ADD_LIBRARY(vtkshp SHP_SRCS)

INSTALL_TARGETS(/lib/vtk vtkshp)
-----------------------------------------------------------
 Also I modify file VTK\Utilities\CMakeLists.txt
 (add shapelib to SUBDIRS before Doxygen directory)
-----------------------------------------------------------
 File VTK/Utilities/CMakeLists.txt:
-----------------------------------------------------------

# build zlib, png, jpeg and shapelib
SUBDIRS(zlib png jpeg shapelib Doxygen)
-----------------------------------------------------------
 Is this simple way correct?
 ---------------------------------------------------------------------------
---------------------------------------------
Build is ok, now I must write some classes for this library.
-----------------------------------------------------------
 File VTK/IO/CMakeLists.txt:
-----------------------------------------------------------
SOURCE_FILES( IO_SRCS
vtkESRI_PolygonZWriter -
...
-----------------------------------------------------------

Now, I need to add the VTK-like C++ wrap for this library. For example,
PolygonZWriter case:
-----------------------------------------------------------
vtkESRI_PolygonZWriter .h file:
-----------------------------------------------------------
class VTK_IO_EXPORT vtkESRI_PolygonZWriter : public vtkPolyDataWriter
{
public:
  static vtkESRI_PolygonZWriter *New();
  vtkTypeMacro(vtkESRI_PolygonZWriter,vtkPolyDataWriter);

protected:
  vtkESRI_PolygonZWriter();
  ~vtkESRI_PolygonZWriter() {};

  void WriteData();

private:
  vtkESRI_PolygonZWriter(const vtkESRI_PolygonZWriter&);  // Not
implemented.
  void operator=(const vtkESRI_PolygonZWriter&);  // Not implemented.
};
-----------------------------------------------------------
.cpp file, general function:
-----------------------------------------------------------
void vtkESRI_PolygonZWriter::WriteData()
{
 vtkPoints *pts;
 vtkCellArray *polys;

  vtkPolyData *input = this->GetInput();
  int numCells = input->GetNumberOfCells();
  vtkCellData *cd= input->GetCellData();
  vtkDataArray *scalars = cd->GetScalars();

  int dataType;
  int numComp;
  if ( scalars && scalars->GetNumberOfTuples() > 0 )
  {
   dataType = scalars->GetDataType();
   numComp = scalars->GetNumberOfComponents();
  }

  polys = input->GetPolys();
  pts = input->GetPoints();
  if (pts == NULL || polys == NULL )
    {
    vtkErrorMacro(<<"No data to write!");
    return;
    }

  if ( this->FileName == NULL)
    {
    vtkErrorMacro(<< "Please specify FileName to write");
    return;
    }

  FILE *fp;
  float n[3],  *v;
  double *x, *y, *z;
  vtkIdType npts;
  vtkIdType *indx;

  if ((fp = fopen(this->FileName, "w")) == NULL)
    {
    vtkErrorMacro(<< "Couldn't open file: " << this->FileName);
    return;
    }
   vtkDebugMacro("Writing ESRI POLYGONZ shape file ");

 SHPHandle hSHPHandle = SHPCreate( this->FileName, SHPT_POLYGONZ );
 DBFHandle psDBF = DBFCreate( this->FileName );
 DBFAddField(psDBF, "Id", FTInteger, 10, 0);
  if ( scalars && scalars->GetNumberOfTuples() > 0 )
   {

  switch (dataType)
  {
   case VTK_BIT:
    DBFAddField(psDBF, scalars->GetName(), FTInteger, 1, 0);
   case VTK_CHAR:
   case VTK_UNSIGNED_CHAR:
    DBFAddField(psDBF, scalars->GetName(), FTInteger, 4, 0);
   break;
   case VTK_SHORT:
   case VTK_UNSIGNED_SHORT:
    DBFAddField(psDBF, scalars->GetName(), FTInteger, 8, 0);
   break;

   case VTK_INT:
   case VTK_ID_TYPE:// currently writing vtkIdType as int.
   case VTK_UNSIGNED_INT:
    DBFAddField(psDBF, scalars->GetName(), FTInteger, 16, 0);
   break;

   case VTK_LONG:
   case VTK_UNSIGNED_LONG:
    DBFAddField(psDBF, scalars->GetName(), FTInteger, 32, 0);
   break;

   case VTK_FLOAT:
    DBFAddField(psDBF, scalars->GetName(), FTDouble, 10, 2);
   break;
   case VTK_DOUBLE:
    DBFAddField(psDBF, scalars->GetName(), FTDouble, 16, 6);
   break;
  }
  }
 int j=0;
 for (polys->InitTraversal(); polys->GetNextCell(npts,indx); )
    {
  // write geometry
  x = new double[npts+1];
  y = new double[npts+1];
  z = new double[npts+1];
  for(int i=0; i<npts; i++)
  {
   v = pts->GetPoint(indx[i]);
   x[i] = v[0];
   y[i] = v[1];
   z[i] = v[2];
  }
  x[npts] = x[0];
  y[npts] = y[0];
  z[npts] = z[0];
  SHPObject *pObject = SHPCreateSimpleObject( SHPT_POLYGONZ, npts+1, x,y,z);
  SHPWriteObject(hSHPHandle, -1,  pObject );
  DBFWriteIntegerAttribute( psDBF, j, 0, j+1);
  SHPDestroyObject(pObject);
  delete x; delete y; delete z;
  // write attribute data
  if ( scalars && scalars->GetNumberOfTuples() > 0 )
  {

   float fv = ((vtkFloatArray *)scalars)->GetValue(j);
   //DBFWriteIntegerAttribute( psDBF, j, 0, j);
   DBFWriteDoubleAttribute( psDBF, j, 1, fv);
   j++;
  }

 }
 SHPClose(hSHPHandle);
 DBFClose(psDBF);
}