QGIS/src/analysis/vector/qgstransectsample.cpp

#include "qgstransectsample.h"
#include "qgsdistancearea.h"
#include "qgsgeometry.h"
#include "qgsspatialindex.h"
#include "qgsvectorfilewriter.h"
#include "qgsvectorlayer.h"
#include <QProgressDialog>
#include <QFileInfo>
#include <stdint.h>
#include "mersenne-twister.h"
#include <limits>

QgsTransectSample::QgsTransectSample( QgsVectorLayer* strataLayer, QString strataIdAttribute, QString minDistanceAttribute, QString nPointsAttribute, DistanceUnits minDistUnits,
                                      QgsVectorLayer* baselineLayer, bool shareBaseline, QString baselineStrataId, const QString& outputPointLayer,
                                      const QString& outputLineLayer, const QString& usedBaselineLayer, double minTransectLength ): mStrataLayer( strataLayer ),
    mStrataIdAttribute( strataIdAttribute ), mMinDistanceAttribute( minDistanceAttribute ), mNPointsAttribute( nPointsAttribute ), mBaselineLayer( baselineLayer ), mShareBaseline( shareBaseline ),
    mBaselineStrataId( baselineStrataId ), mOutputPointLayer( outputPointLayer ), mOutputLineLayer( outputLineLayer ), mUsedBaselineLayer( usedBaselineLayer ),
    mMinDistanceUnits( minDistUnits ), mMinTransectLength( minTransectLength )
{
}

QgsTransectSample::QgsTransectSample()
{
}

QgsTransectSample::~QgsTransectSample()
{
}

int QgsTransectSample::createSample( QProgressDialog* pd )
{
  Q_UNUSED( pd );

  if ( !mStrataLayer || !mStrataLayer->isValid() )
  {
    return 1;
  }

  if ( !mBaselineLayer || !mBaselineLayer->isValid() )
  {
    return 2;
  }

  //stratum id is not necessarily an integer
  QVariant::Type stratumIdType = QVariant::Int;
  if ( !mStrataIdAttribute.isEmpty() )
  {
    stratumIdType = mStrataLayer->pendingFields().field( mStrataIdAttribute ).type();
  }

  //create vector file writers for output
  QgsFields outputPointFields;
  outputPointFields.append( QgsField( "id", stratumIdType ) );
  outputPointFields.append( QgsField( "station_id", QVariant::Int ) );
  outputPointFields.append( QgsField( "stratum_id", stratumIdType ) );
  outputPointFields.append( QgsField( "station_code", QVariant::String ) );
  outputPointFields.append( QgsField( "start_lat", QVariant::Double ) );
  outputPointFields.append( QgsField( "start_long", QVariant::Double ) );

  QgsVectorFileWriter outputPointWriter( mOutputPointLayer, "utf-8", outputPointFields, QGis::WKBPoint,
                                         &( mStrataLayer->crs() ) );
  if ( outputPointWriter.hasError() != QgsVectorFileWriter::NoError )
  {
    return 3;
  }

  outputPointFields.append( QgsField( "bearing", QVariant::Double ) ); //add bearing attribute for lines
  QgsVectorFileWriter outputLineWriter( mOutputLineLayer, "utf-8", outputPointFields, QGis::WKBLineString,
                                        &( mStrataLayer->crs() ) );
  if ( outputLineWriter.hasError() != QgsVectorFileWriter::NoError )
  {
    return 4;
  }

  QgsFields usedBaselineFields;
  usedBaselineFields.append( QgsField( "stratum_id", stratumIdType ) );
  usedBaselineFields.append( QgsField( "ok", QVariant::String ) );
  QgsVectorFileWriter usedBaselineWriter( mUsedBaselineLayer, "utf-8", usedBaselineFields, QGis::WKBLineString,
                                          &( mStrataLayer->crs() ) );
  if ( usedBaselineWriter.hasError() != QgsVectorFileWriter::NoError )
  {
    return 5;
  }

  //debug: write clipped buffer bounds with stratum id to same directory as out_point
  QFileInfo outputPointInfo( mOutputPointLayer );
  QString bufferClipLineOutput = outputPointInfo.absolutePath() + "/out_buffer_clip_line.shp";
  QgsFields bufferClipLineFields;
  bufferClipLineFields.append( QgsField( "id", stratumIdType ) );
  QgsVectorFileWriter bufferClipLineWriter( bufferClipLineOutput, "utf-8", bufferClipLineFields, QGis::WKBLineString, &( mStrataLayer->crs() ) );

  //configure distanceArea depending on minDistance units and output CRS
  QgsDistanceArea distanceArea;
  distanceArea.setSourceCrs( mStrataLayer->crs().srsid() );
  if ( mMinDistanceUnits == Meters )
  {
    distanceArea.setEllipsoidalMode( true );
  }
  else
  {
    distanceArea.setEllipsoidalMode( false );
  }

  //possibility to transform output points to lat/long
  QgsCoordinateTransform toLatLongTransform( mStrataLayer->crs(), QgsCoordinateReferenceSystem( 4326, QgsCoordinateReferenceSystem::EpsgCrsId ) );

  //init random number generator
  mt_srand( QTime::currentTime().msec() );

  QgsFeatureRequest fr;
  fr.setSubsetOfAttributes( QStringList() << mStrataIdAttribute << mMinDistanceAttribute << mNPointsAttribute, mStrataLayer->pendingFields() );
  QgsFeatureIterator strataIt = mStrataLayer->getFeatures( fr );

  QgsFeature fet;
  int nTotalTransects = 0;
  int nFeatures = 0;

  if ( pd )
  {
    pd->setMaximum( mStrataLayer->featureCount() );
  }

  while ( strataIt.nextFeature( fet ) )
  {
    if ( pd )
    {
      pd->setValue( nFeatures );
    }
    if ( pd && pd->wasCanceled() )
    {
      break;
    }


    QgsGeometry* strataGeom = fet.geometry();
    if ( !strataGeom )
    {
      continue;
    }

    //find baseline for strata
    bool strataIdOk = true;
    QVariant strataId = fet.attribute( mStrataIdAttribute );
    QgsGeometry* baselineGeom = findBaselineGeometry( strataIdOk ? strataId : -1 );
    if ( !baselineGeom )
    {
      continue;
    }

    double minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
    double minDistanceLayerUnits = minDistance;
    //if minDistance is in meters and the data in degrees, we need to apply a rough conversion for the buffer distance
    double bufferDist = minDistance;
    if ( mMinDistanceUnits == Meters && mStrataLayer->crs().mapUnits() == QGis::DecimalDegrees )
    {
      bufferDist = minDistance / 111319.9;
      minDistanceLayerUnits = bufferDist;
    }

    QgsGeometry* clippedBaseline = strataGeom->intersection( baselineGeom );
    if ( !clippedBaseline || clippedBaseline->wkbType() == QGis::WKBUnknown )
    {
      delete clippedBaseline;
      continue;
    }
    QgsGeometry* bufferLineClipped = clipBufferLine( strataGeom, clippedBaseline, bufferDist );
    if ( !bufferLineClipped )
    {
      delete clippedBaseline;
      continue;
    }

    //save clipped baseline to file
    QgsFeature blFeature;
    blFeature.setGeometry( *clippedBaseline );
    blFeature.setAttribute( "stratum_id", strataId );
    blFeature.setAttribute( "ok", "f" );
    usedBaselineWriter.addFeature( blFeature );

    //start loop to create random points along the baseline
    int nTransects = fet.attribute( mNPointsAttribute ).toInt();
    int nCreatedTransects = 0;
    int nIterations = 0;
    int nMaxIterations = nTransects * 50;

    QgsSpatialIndex sIndex; //to check minimum distance
    QMap< QgsFeatureId, QgsGeometry* > lineFeatureMap;

    while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
    {
      double randomPosition = (( double )mt_rand() / RAND_MAX ) * clippedBaseline->length();
      QgsGeometry* samplePoint = clippedBaseline->interpolate( randomPosition );
      ++nIterations;
      if ( !samplePoint )
      {
        continue;
      }
      QgsPoint sampleQgsPoint = samplePoint->asPoint();
      QgsPoint latLongSamplePoint = toLatLongTransform.transform( sampleQgsPoint );

      QgsFeature samplePointFeature;
      samplePointFeature.setGeometry( samplePoint );
      samplePointFeature.setAttribute( "id", nTotalTransects + 1 );
      samplePointFeature.setAttribute( "station_id", nCreatedTransects + 1 );
      samplePointFeature.setAttribute( "stratum_id", strataId );
      samplePointFeature.setAttribute( "station_code", strataId.toString() + "_" + QString::number( nCreatedTransects + 1 ) );
      samplePointFeature.setAttribute( "start_lat", latLongSamplePoint.y() );
      samplePointFeature.setAttribute( "start_long", latLongSamplePoint.x() );

      //find closest point on clipped buffer line
      QgsPoint minDistPoint;

      int afterVertex;
      if ( bufferLineClipped->closestSegmentWithContext( sampleQgsPoint, minDistPoint, afterVertex ) < 0 )
      {
        continue;
      }

      //bearing between sample point and min dist point (transect direction)
      double bearing = distanceArea.bearing( sampleQgsPoint, minDistPoint ) / M_PI * 180.0;

      QgsPolyline sampleLinePolyline;
      QgsPoint ptFarAway( sampleQgsPoint.x() + ( minDistPoint.x() - sampleQgsPoint.x() ) * 1000000,
                          sampleQgsPoint.y() + ( minDistPoint.y() - sampleQgsPoint.y() ) * 1000000 );
      QgsPolyline lineFarAway;
      lineFarAway << sampleQgsPoint << ptFarAway;
      QgsGeometry* lineFarAwayGeom = QgsGeometry::fromPolyline( lineFarAway );
      QgsGeometry* lineClipStratum = lineFarAwayGeom->intersection( strataGeom );
      if ( !lineClipStratum )
      {
        delete lineFarAwayGeom; delete lineClipStratum;
        continue;
      }

      //cancel if distance between sample point and line is too large (line does not start at point
      if ( lineClipStratum->distance( *samplePoint ) > 0.000001 )
      {
        delete lineFarAwayGeom; delete lineClipStratum;
        continue;
      }

      //if lineClipStratum is a multiline, take the part line closest to sampleQgsPoint
      if ( lineClipStratum->wkbType() == QGis::WKBMultiLineString
           || lineClipStratum->wkbType() == QGis::WKBMultiLineString25D )
      {
        QgsGeometry* singleLine = closestMultilineElement( sampleQgsPoint, lineClipStratum );
        if ( singleLine )
        {
          delete lineClipStratum;
          lineClipStratum = singleLine;
        }
      }

      //cancel if length of lineClipStratum is too small
      double transectLength = distanceArea.measure( lineClipStratum );
      if ( transectLength < mMinTransectLength )
      {
        delete lineFarAwayGeom; delete lineClipStratum;
        continue;
      }

      //search closest existing profile. Cancel if dist < minDist
      if ( otherTransectWithinDistance( lineClipStratum, minDistanceLayerUnits, minDistance, sIndex, lineFeatureMap, distanceArea ) )
      {
        delete lineFarAwayGeom; delete lineClipStratum;
        continue;
      }

      QgsFeatureId fid( nCreatedTransects );
      QgsFeature sampleLineFeature( fid );
      sampleLineFeature.setGeometry( lineClipStratum );
      sampleLineFeature.setAttribute( "id", nTotalTransects + 1 );
      sampleLineFeature.setAttribute( "station_id", nCreatedTransects + 1 );
      sampleLineFeature.setAttribute( "stratum_id", strataId );
      sampleLineFeature.setAttribute( "station_code", strataId.toString() + "_" + QString::number( nCreatedTransects + 1 ) );
      sampleLineFeature.setAttribute( "start_lat", latLongSamplePoint.y() );
      sampleLineFeature.setAttribute( "start_long", latLongSamplePoint.x() );
      sampleLineFeature.setAttribute( "bearing", bearing );
      outputLineWriter.addFeature( sampleLineFeature );

      //add point to file writer here.
      //It can only be written if the corresponding transect has been as well
      outputPointWriter.addFeature( samplePointFeature );

      sIndex.insertFeature( sampleLineFeature );
      lineFeatureMap.insert( fid, sampleLineFeature.geometryAndOwnership() );

      delete lineFarAwayGeom;
      ++nTotalTransects;
      ++nCreatedTransects;
    }
    delete clippedBaseline;

    QgsFeature bufferClipFeature;
    bufferClipFeature.setGeometry( bufferLineClipped );
    bufferClipFeature.setAttribute( "id" , strataId );
    bufferClipLineWriter.addFeature( bufferClipFeature );
    //delete bufferLineClipped;

    //delete all line geometries in spatial index
    QMap< QgsFeatureId, QgsGeometry* >::iterator featureMapIt = lineFeatureMap.begin();
    for ( ; featureMapIt != lineFeatureMap.end(); ++featureMapIt )
    {
      delete( featureMapIt.value() );
    }
    lineFeatureMap.clear();
    delete baselineGeom;

    ++nFeatures;
  }

  if ( pd )
  {
    pd->setValue( mStrataLayer->featureCount() );
  }

  return 0;
}

QgsGeometry* QgsTransectSample::findBaselineGeometry( QVariant strataId )
{
  if ( !mBaselineLayer )
  {
    return 0;
  }

  QgsFeatureIterator baseLineIt = mBaselineLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QStringList( mBaselineStrataId ), mBaselineLayer->pendingFields() ) );
  QgsFeature fet;

  while ( baseLineIt.nextFeature( fet ) ) //todo: cache this in case there are many baslines
  {
    if ( strataId == fet.attribute( mBaselineStrataId ) || mShareBaseline )
    {
      return fet.geometryAndOwnership();
    }
  }
  return 0;
}

bool QgsTransectSample::otherTransectWithinDistance( QgsGeometry* geom, double minDistLayerUnit, double minDistance, QgsSpatialIndex& sIndex,
    const QMap< QgsFeatureId, QgsGeometry* >& lineFeatureMap, QgsDistanceArea& da )
{
  if ( !geom )
  {
    return false;
  }

  QgsGeometry* buffer = geom->buffer( minDistLayerUnit, 8 );
  if ( !buffer )
  {
    return false;
  }
  QgsRectangle rect = buffer->boundingBox();
  QList<QgsFeatureId> lineIdList = sIndex.intersects( rect );

  QList<QgsFeatureId>::const_iterator lineIdIt = lineIdList.constBegin();
  for ( ; lineIdIt != lineIdList.constEnd(); ++lineIdIt )
  {
    const QMap< QgsFeatureId, QgsGeometry* >::const_iterator idMapIt = lineFeatureMap.find( *lineIdIt );
    if ( idMapIt != lineFeatureMap.constEnd() )
    {
      double dist = 0;
      QgsPoint pt1, pt2;
      closestSegmentPoints( *geom, *( idMapIt.value() ), dist, pt1, pt2 );
      dist = da.measureLine( pt1, pt2 ); //convert degrees to meters if necessary

      if ( dist < minDistance )
      {
        delete buffer;
        return true;
      }
    }
  }

  delete buffer;
  return false;
}

bool QgsTransectSample::closestSegmentPoints( QgsGeometry& g1, QgsGeometry& g2, double& dist, QgsPoint& pt1, QgsPoint& pt2 )
{
  QGis::WkbType t1 = g1.wkbType();
  if ( t1 != QGis::WKBLineString && t1 != QGis::WKBLineString25D )
  {
    return false;
  }

  QGis::WkbType t2 = g2.wkbType();
  if ( t2 != QGis::WKBLineString && t2 != QGis::WKBLineString25D )
  {
    return false;
  }

  QgsPolyline pl1 = g1.asPolyline();
  QgsPolyline pl2 = g2.asPolyline();

  if ( pl1.size() < 2 || pl2.size() < 2 )
  {
    return false;
  }

  QgsPoint p11 = pl1.at( 0 );
  QgsPoint p12 = pl1.at( 1 );
  QgsPoint p21 = pl2.at( 0 );
  QgsPoint p22 = pl2.at( 1 );

  double p1x = p11.x();
  double p1y = p11.y();
  double v1x = p12.x() - p11.x();
  double v1y = p12.y() - p11.y();
  double p2x = p21.x();
  double p2y = p21.y();
  double v2x = p22.x() - p21.x();
  double v2y = p22.y() - p21.y();

  double denominatorU = v2x * v1y - v2y * v1x;
  double denominatorT = v1x * v2y - v1y * v2x;

  if ( qgsDoubleNear( denominatorU, 0 ) || qgsDoubleNear( denominatorT, 0 ) )
  {
    //lines are parallel
    //project all points on the other segment and take the one with the smallest distance
    QgsPoint minDistPoint1;
    double d1 = p11.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint1 );
    QgsPoint minDistPoint2;
    double d2 = p12.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), minDistPoint2 );
    QgsPoint minDistPoint3;
    double d3 = p21.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint3 );
    QgsPoint minDistPoint4;
    double d4 = p22.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), minDistPoint4 );

    if ( d1 <= d2 && d1 <= d3 && d1 <= d4 )
    {
      dist = sqrt( d1 ); pt1 = p11; pt2 = minDistPoint1;
      return true;
    }
    else if ( d2 <= d1 && d2 <= d3 && d2 <= d4 )
    {
      dist = sqrt( d2 );  pt1 = p12; pt2 = minDistPoint2;
      return true;
    }
    else if ( d3 <= d1 && d3 <= d2 && d3 <= d4 )
    {
      dist = sqrt( d3 ); pt1 = p21; pt2 = minDistPoint3;
      return true;
    }
    else
    {
      dist = sqrt( d4 ); pt1 = p21; pt2 = minDistPoint4;
      return true;
    }
  }

  double u = ( p1x * v1y - p1y * v1x - p2x * v1y + p2y * v1x ) / denominatorU;
  double t = ( p2x * v2y - p2y * v2x - p1x * v2y + p1y * v2x ) / denominatorT;

  if ( u >= 0 && u <= 1.0 && t >= 0 && t <= 1.0 )
  {
    dist = 0;
    pt1.setX( p2x + u * v2x );
    pt1.setY( p2y + u * v2y );
    pt2 = pt1;
    dist = 0;
    return true;
  }

  if ( t > 1.0 )
  {
    pt1.setX( p12.x() );
    pt1.setY( p12.y() );
  }
  else if ( t < 0.0 )
  {
    pt1.setX( p11.x() );
    pt1.setY( p11.y() );
  }
  if ( u > 1.0 )
  {
    pt2.setX( p22.x() );
    pt2.setY( p22.y() );
  }
  if ( u < 0.0 )
  {
    pt2.setX( p21.x() );
    pt2.setY( p21.y() );
  }
  if ( t >= 0.0 && t <= 1.0 )
  {
    //project pt2 onto g1
    pt2.sqrDistToSegment( p11.x(), p11.y(), p12.x(), p12.y(), pt1 );
  }
  if ( u >= 0.0 && u <= 1.0 )
  {
    //project pt1 onto g2
    pt1.sqrDistToSegment( p21.x(), p21.y(), p22.x(), p22.y(), pt2 );
  }

  dist = sqrt( pt1.sqrDist( pt2 ) );
  return true;
}

QgsGeometry* QgsTransectSample::closestMultilineElement( const QgsPoint& pt, QgsGeometry* multiLine )
{
  if ( !multiLine || ( multiLine->wkbType() != QGis::WKBMultiLineString
                       && multiLine->wkbType() != QGis::WKBMultiLineString25D ) )
  {
    return 0;
  }

  double minDist = DBL_MAX;
  double currentDist = 0;
  QgsGeometry* currentLine = 0;
  QgsGeometry* closestLine = 0;
  QgsGeometry* pointGeom = QgsGeometry::fromPoint( pt );

  QgsMultiPolyline multiPolyline = multiLine->asMultiPolyline();
  QgsMultiPolyline::const_iterator it = multiPolyline.constBegin();
  for ( ; it != multiPolyline.constEnd(); ++it )
  {
    currentLine = QgsGeometry::fromPolyline( *it );
    currentDist = pointGeom->distance( *currentLine );
    if ( currentDist < minDist )
    {
      minDist = currentDist;
      closestLine = currentLine;
    }
    else
    {
      delete currentLine;
    }
  }

  delete pointGeom;
  return closestLine;
}

QgsGeometry* QgsTransectSample::clipBufferLine( QgsGeometry* stratumGeom, QgsGeometry* clippedBaseline, double tolerance )
{
  if ( !stratumGeom || !clippedBaseline || clippedBaseline->wkbType() == QGis::WKBUnknown )
  {
    return 0;
  }

  double currentBufferDist = tolerance;
  int maxLoops = 10;

  for ( int i = 0; i < maxLoops; ++i )
  {
    //loop with tolerance: create buffer, convert buffer to line, clip line by stratum, test if result is (single) line
    QgsGeometry* clipBaselineBuffer = clippedBaseline->buffer( currentBufferDist, 8 );
    if ( !clipBaselineBuffer )
    {
      delete clipBaselineBuffer;
      continue;
    }

    //it is also possible that clipBaselineBuffer is a multipolygon
    QgsGeometry* bufferLine = 0; //buffer line or multiline
    QgsGeometry* bufferLineClipped = 0;
    QgsMultiPolyline mpl;
    if ( clipBaselineBuffer->isMultipart() )
    {
      QgsMultiPolygon bufferMultiPolygon = clipBaselineBuffer->asMultiPolygon();
      if ( bufferMultiPolygon.size() < 1 )
      {
        delete clipBaselineBuffer;
        continue;
      }

      for ( int j = 0; j < bufferMultiPolygon.size(); ++j )
      {
        int size = bufferMultiPolygon.at( j ).size();
        for ( int k = 0; k < size; ++k )
        {
          mpl.append( bufferMultiPolygon.at( j ).at( k ) );
        }
      }
      bufferLine = QgsGeometry::fromMultiPolyline( mpl );
    }
    else
    {
      QgsPolygon bufferPolygon = clipBaselineBuffer->asPolygon();
      if ( bufferPolygon.size() < 1 )
      {
        delete clipBaselineBuffer;
        continue;
      }

      int size = bufferPolygon.size();
      for ( int j = 0; j < size; ++j )
      {
        mpl.append( bufferPolygon[j] );
      }
      bufferLine = QgsGeometry::fromMultiPolyline( mpl );
    }
    bufferLineClipped = bufferLine->intersection( stratumGeom );

    if ( bufferLineClipped && bufferLineClipped->type() == QGis::Line )
    {
      //if stratumGeom is a multipolygon, bufferLineClipped must intersect each part
      bool bufferLineClippedIntersectsStratum = true;
      if ( stratumGeom->wkbType() == QGis::WKBMultiPolygon || stratumGeom->wkbType() == QGis::WKBMultiPolygon25D )
      {
        QVector<QgsPolygon> multiPoly = stratumGeom->asMultiPolygon();
        QVector<QgsPolygon>::const_iterator multiIt = multiPoly.constBegin();
        for ( ; multiIt != multiPoly.constEnd(); ++multiIt )
        {
          QgsGeometry* poly = QgsGeometry::fromPolygon( *multiIt );
          if ( !poly->intersects( bufferLineClipped ) )
          {
            bufferLineClippedIntersectsStratum = false;
            delete poly;
            break;
          }
          delete poly;
        }
      }

      if ( bufferLineClippedIntersectsStratum )
      {
        return bufferLineClipped;
      }
    }

    delete bufferLineClipped; delete clipBaselineBuffer; delete bufferLine;
    currentBufferDist /= 2;
  }

  return 0; //no solution found even with reduced tolerances
}