qgsgrassprovider.cpp

一个非常好的GIS开源新版本

/***************************************************************************
    qgsgrassprovider.cpp -  Data provider for GRASS format
                             -------------------
    begin                : March, 2004
    copyright            : (C) 2004 by Gary E.Sherman, Radim Blazek
    email                : [e-mail], [e-mail]
 ***************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
/* $Id: qgsgrassprovider.cpp 8350 2008-04-14 21:02:22Z jef $ */

#include <string.h>
#include <iostream>
#include <vector>
#include <cfloat>

#include <Q3CString>
#include <QDir>
#include <QString>
#include <QDateTime>
#include <QMessageBox>
#include <QTextCodec>

#include "qgis.h"
#include "qgsapplication.h"
#include "qgsdataprovider.h"
#include "qgsfeature.h"
#include "qgsfield.h"
#include "qgslogger.h"
#include "qgsrect.h"
#include "qgsspatialrefsys.h"

#ifdef _MSC_VER
// enable grass prototypes
#define __STDC__ 1
#endif

extern "C" {
#include <grass/gprojects.h>
#include <grass/gis.h>
#include <grass/dbmi.h>
#include <grass/Vect.h>
}

#ifdef _MSC_VER
#undef __STDC__
#endif

#include "qgsgrass.h"
#include "qgsgrassprovider.h"

std::vector<GLAYER> QgsGrassProvider::mLayers;
std::vector<GMAP> QgsGrassProvider::mMaps;


static QString GRASS_KEY = "grass"; // XXX verify this
static QString GRASS_DESCRIPTION = "Grass provider"; // XXX verify this



QgsGrassProvider::QgsGrassProvider(QString uri)
    : QgsVectorDataProvider(uri)
{
#ifdef QGISDEBUG
  std::cerr << "QgsGrassProvider URI: " << uri.toLocal8Bit().data() << std::endl;
#endif

  QgsGrass::init();

  QTime time;
  time.start();

  mValid = false;

  // Parse URI 
  QDir dir ( uri );  // it is not a directory in fact
  QString myURI = dir.path();  // no dupl '/'

  mLayer = dir.dirName();
  myURI = myURI.left( dir.path().findRev('/') );
  dir = QDir(myURI);
  mMapName = dir.dirName();
  dir.cdUp(); 
  mMapset = dir.dirName();
  dir.cdUp(); 
  mLocation = dir.dirName();
  dir.cdUp(); 
  mGisdbase = dir.path();

#ifdef QGISDEBUG
  std::cerr << "gisdbase: " << mGisdbase.toLocal8Bit().data() << std::endl;
  std::cerr << "location: " << mLocation.toLocal8Bit().data() << std::endl;
  std::cerr << "mapset: "   << mMapset.toLocal8Bit().data() << std::endl;
  std::cerr << "mapName: "  << mMapName.toLocal8Bit().data() << std::endl;
  std::cerr << "layer: "    << mLayer.toLocal8Bit().data() << std::endl;
#endif

  /* Parse Layer, supported layers <field>_point, <field>_line, <field>_area
  *  Layer is opened even if it is empty (has no features) 
  */
  mLayerField = -1;       
  if ( mLayer.compare("boundary") == 0 ) { // currently not used
    mLayerType = BOUNDARY;
    mGrassType = GV_BOUNDARY;
  } else if ( mLayer.compare("centroid") == 0 ) { // currently not used
    mLayerType = CENTROID;
    mGrassType = GV_CENTROID;
  } else {
    mLayerField = grassLayer ( mLayer );
    if ( mLayerField == -1 ) {
      std::cerr << "Invalid layer name, no underscore found: " << mLayer.toLocal8Bit().data() << std::endl;
      return;
    }

    mGrassType = grassLayerType ( mLayer );

    if ( mGrassType == GV_POINT ) {
      mLayerType = POINT;
    } else if ( mGrassType == GV_LINES ) {
      mLayerType = LINE;
    } else if ( mGrassType == GV_AREA ) {
      mLayerType = POLYGON;
    } else {
      std::cerr << "Invalid layer name, wrong type: " << mLayer.toLocal8Bit().data() << std::endl;
      return;
    }

  }
#ifdef QGISDEBUG
  std::cerr << "mLayerField: " << mLayerField << std::endl;
  std::cerr << "mLayerType: " << mLayerType << std::endl;
#endif

  if ( mLayerType == BOUNDARY || mLayerType == CENTROID ) {
    std::cerr << "Layer type not supported." << std::endl;
    return;
  }

  // Set QGIS type
  switch ( mLayerType ) {
  case POINT:
  case CENTROID:
    mQgisType = QGis::WKBPoint; 
    break;
  case LINE:
  case BOUNDARY:
    mQgisType = QGis::WKBLineString; 
    break;
  case POLYGON:
    mQgisType = QGis::WKBPolygon;
    break;
  }

  mLayerId = openLayer(mGisdbase, mLocation, mMapset, mMapName, mLayerField);
  if ( mLayerId < 0 ) {
    std::cerr << "Cannot open GRASS layer:" << myURI.toLocal8Bit().data() << std::endl;
    return;
  }
#ifdef QGISDEBUG
  std::cerr << "mLayerId: " << mLayerId << std::endl;
#endif

  mMap = layerMap(mLayerId);

  // Getting the total number of features in the layer
  mNumberFeatures = 0;
  mCidxFieldIndex = -1;
  if ( mLayerField >= 0 ) {
    mCidxFieldIndex = Vect_cidx_get_field_index ( mMap, mLayerField);
    if ( mCidxFieldIndex >= 0 ) {
      mNumberFeatures = Vect_cidx_get_type_count ( mMap, mLayerField, mGrassType );
      mCidxFieldNumCats = Vect_cidx_get_num_cats_by_index ( mMap, mCidxFieldIndex );
    }
  } else {
    // TODO nofield layers
    mNumberFeatures = 0;
    mCidxFieldNumCats = 0;
  }
  mNextCidx = 0;

#ifdef QGISDEBUG
  std::cerr << "mNumberFeatures = " << mNumberFeatures << " mCidxFieldIndex = " << mCidxFieldIndex
    << " mCidxFieldNumCats = " << mCidxFieldNumCats << std::endl;
#endif


  // Create selection array
  mSelectionSize = allocateSelection ( mMap, &mSelection );
  resetSelection(1); // TODO ? - where what reset

  mMapVersion = mMaps[mLayers[mLayerId].mapId].version;

  // Init structures
  mPoints = Vect_new_line_struct ();
  mCats = Vect_new_cats_struct ();
  mList = Vect_new_list ();

  mValid = true;

#ifdef QGISDEBUG
  std::cerr << "New GRASS layer opened, time (ms): " << time.elapsed() << std::endl;
#endif
}

void QgsGrassProvider::update ( void )
{
#ifdef QGISDEBUG
  std::cerr << "*** QgsGrassProvider::update ***" << std::endl;
#endif

  mValid = false;

  if ( ! mMaps[mLayers[mLayerId].mapId].valid ) return;

  // Getting the total number of features in the layer
  // It may happen that the field disappeares from the map (deleted features, new map without that field)
  mNumberFeatures = 0;
  mCidxFieldIndex = -1;
  if ( mLayerField >= 0 ) {
    mCidxFieldIndex = Vect_cidx_get_field_index ( mMap, mLayerField);
    if ( mCidxFieldIndex >= 0 ) {
      mNumberFeatures = Vect_cidx_get_type_count ( mMap, mLayerField, mGrassType );
      mCidxFieldNumCats = Vect_cidx_get_num_cats_by_index ( mMap, mCidxFieldIndex );
    }
  } else {
    // TODO nofield layers
    mNumberFeatures = 0;
    mCidxFieldNumCats = 0;
  }
  mNextCidx = 0;

#ifdef QGISDEBUG
  std::cerr << "mNumberFeatures = " << mNumberFeatures << " mCidxFieldIndex = " << mCidxFieldIndex
    << " mCidxFieldNumCats = " << mCidxFieldNumCats << std::endl;
#endif

  // Create selection array
  if ( mSelection ) free ( mSelection );
  mSelectionSize = allocateSelection ( mMap, &mSelection );
  resetSelection(1); 

  mMapVersion = mMaps[mLayers[mLayerId].mapId].version;

  mValid = true;
}

int QgsGrassProvider::allocateSelection( struct Map_info *map, char **selection )
{
  int size;
#ifdef QGISDEBUG
  std::cerr << "QgsGrassProvider::allocateSellection" << std::endl;
#endif

  int nlines = Vect_get_num_lines ( map );
  int nareas = Vect_get_num_areas ( map );

  if ( nlines > nareas ) {
    size = nlines + 1;
  } else {
    size = nareas + 1;
  }
#ifdef QGISDEBUG
  std::cerr << "nlines = " << nlines << " nareas = " << nareas << " size = " << size << std::endl;
#endif

  *selection = (char *) malloc ( size );

  return size;
}

QgsGrassProvider::~QgsGrassProvider()
{
#ifdef QGISDEBUG
  std::cerr << "QgsGrassProvider::~QgsGrassProvider()" << std::endl;
#endif
  closeLayer ( mLayerId );
}


QString QgsGrassProvider::storageType() const
{
  return "GRASS (Geographic Resources Analysis and Support System) file";
}

bool QgsGrassProvider::getNextFeature(QgsFeature& feature)
{
  int cat, type, id;
  unsigned char *wkb;
  int wkbsize;

#if QGISDEBUG > 3
  std::cout << "QgsGrassProvider::getNextFeature()" << std::endl;
#endif

  if ( isEdited() || isFrozen() || !mValid )
    return false;

  if ( mCidxFieldIndex < 0 ) return false; // No features, no features in this layer

  // Get next line/area id
  int found = 0;
  while ( mNextCidx < mCidxFieldNumCats ) {
    Vect_cidx_get_cat_by_index ( mMap, mCidxFieldIndex, mNextCidx++, &cat, &type, &id );
    // Warning: selection array is only of type line/area of current layer -> check type first

    if ( !(type & mGrassType) ) continue;
    if ( !mSelection[id] ) continue;
    found = 1;
    break;
  }
  if ( !found ) return false; // No more features
#if QGISDEBUG > 3
  std::cout << "cat = " << cat << " type = " << type << " id = " << id << std::endl;
#endif

  feature = QgsFeature(id);

  // TODO int may be 64 bits (memcpy)
  if ( type & (GV_POINTS | GV_LINES) ) { /* points or lines */
    Vect_read_line ( mMap, mPoints, mCats, id);
    int npoints = mPoints->n_points;

    if ( type & GV_POINTS ) {
      wkbsize = 1 + 4 + 2*8;
    } else { // GV_LINES
      wkbsize = 1+4+4+npoints*2*8;
    }	    
    wkb = new unsigned char[wkbsize];
    unsigned char *wkbp = wkb;
    wkbp[0] = (unsigned char) QgsApplication::endian();
    wkbp += 1;

    /* WKB type */
    memcpy (wkbp, &mQgisType, 4);
    wkbp += 4;

    /* number of points */
    if ( type & GV_LINES ) {
      memcpy (wkbp, &npoints, 4);
      wkbp += 4;
    }

    for ( int i = 0; i < npoints; i++ ) {
      memcpy (wkbp, &(mPoints->x[i]), 8);
      memcpy (wkbp+8, &(mPoints->y[i]), 8);
      wkbp += 16;
    }
  } else { // GV_AREA
    Vect_get_area_points ( mMap, id, mPoints );
    int npoints = mPoints->n_points;

    wkbsize = 1+4+4+4+npoints*2*8; // size without islands
    wkb = new unsigned char[wkbsize];
    wkb[0] = (unsigned char) QgsApplication::endian();
    int offset = 1;

    /* WKB type */
    memcpy ( wkb+offset, &mQgisType, 4);
    offset += 4;

    /* Number of rings */
    int nisles = Vect_get_area_num_isles ( mMap, id );
    int nrings = 1 + nisles; 
    memcpy (wkb+offset, &nrings, 4);
    offset += 4;

    /* Outer ring */
    memcpy (wkb+offset, &npoints, 4);
    offset += 4;
    for ( int i = 0; i < npoints; i++ ) {
      memcpy (wkb+offset, &(mPoints->x[i]), 8);
      memcpy (wkb+offset+8, &(mPoints->y[i]), 8);
      offset += 16;
    }

    /* Isles */
    for ( int i = 0; i < nisles; i++ ) {
      Vect_get_isle_points ( mMap, Vect_get_area_isle (mMap, id, i), mPoints );
      npoints = mPoints->n_points;

      // add space
      wkbsize += 4+npoints*2*8;
      wkb = (unsigned char *) realloc (wkb, wkbsize);

      memcpy (wkb+offset, &npoints, 4);
      offset += 4;
      for ( int i = 0; i < npoints; i++ ) {
        memcpy (wkb+offset, &(mPoints->x[i]), 8);
        memcpy (wkb+offset+8, &(mPoints->y[i]), 8);
        offset += 16;
      }
    }
  }

  feature.setGeometryAndOwnership(wkb, wkbsize);

  setFeatureAttributes( mLayerId, cat, &feature, mAttributesToFetch);  

  return true;
}

void QgsGrassProvider::resetSelection( bool sel)
{
#ifdef QGISDEBUG
  std::cout << "QgsGrassProvider::resetSelection()" << std::endl;
#endif
  if ( !mValid ) return;
  memset ( mSelection, (int) sel, mSelectionSize );
  mNextCidx = 0;
}

void QgsGrassProvider::select(QgsAttributeList fetchAttributes,
                              QgsRect rect,
                              bool fetchGeometry,
                              bool useIntersect)
{
  mAttributesToFetch = fetchAttributes;
  mFetchGeom = fetchGeometry;

  if ( isEdited() || isFrozen() || !mValid )
    return;

  // check if outdated and update if necessary
  int mapId = mLayers[mLayerId].mapId;
  if ( mapOutdated(mapId) ) 
  {
    updateMap ( mapId );
  }
  if ( mMapVersion < mMaps[mapId].version ) 
  {
    update();
  }
  if ( attributesOutdated(mapId) ) 
  {
    loadAttributes (mLayers[mLayerId]);
  }

  //no selection rectangle - use all features
  if(rect.isEmpty())
  {
    resetSelection(1);
    return;
  }

  //apply selection rectangle
  resetSelection(0);

  if ( !useIntersect ) 
  { // select by bounding boxes only
    BOUND_BOX box;
    box.N = rect.yMax(); box.S = rect.yMin(); 
    box.E = rect.xMax(); box.W = rect.xMin(); 
    box.T = PORT_DOUBLE_MAX; box.B = -PORT_DOUBLE_MAX; 
    if ( mLayerType == POINT || mLayerType == CENTROID || mLayerType == LINE || mLayerType == BOUNDARY ) 
    {
      Vect_select_lines_by_box(mMap, &box, mGrassType, mList);
    } 
    else if ( mLayerType == POLYGON ) 
    {
      Vect_select_areas_by_box(mMap, &box, mList);
    }

  } 
  else 
  { // check intersection
    struct line_pnts *Polygon;

    Polygon = Vect_new_line_struct();

    Vect_append_point( Polygon, rect.xMin(), rect.yMin(), 0);
    Vect_append_point( Polygon, rect.xMax(), rect.yMin(), 0);
    Vect_append_point( Polygon, rect.xMax(), rect.yMax(), 0);