qgsrasterlayerrenderer.cpp

Go to the documentation of this file.

/***************************************************************************
  qgsrasterlayerrenderer.cpp
  --------------------------------------
  Date                 : December 2013
  Copyright            : (C) 2013 by Martin Dobias
  Email                : wonder dot sk at gmail dot com
 ***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgsrasterlayerrenderer.h"
 
#include <memory>
 
#include "qgsapplication.h"
#include "qgselevationmap.h"
#include "qgsexception.h"
#include "qgsgdalutils.h"
#include "qgsinterval.h"
#include "qgsmapclippingutils.h"
#include "qgsmessagelog.h"
#include "qgsrasterdataprovider.h"
#include "qgsrasterdrawer.h"
#include "qgsrasteriterator.h"
#include "qgsrasterlabeling.h"
#include "qgsrasterlayer.h"
#include "qgsrasterlayerelevationproperties.h"
#include "qgsrasterlayertemporalproperties.h"
#include "qgsrasterlayerutils.h"
#include "qgsrasternuller.h"
#include "qgsrasterpipe.h"
#include "qgsrasterprojector.h"
#include "qgsrasterrenderer.h"
#include "qgsrasterresamplefilter.h"
#include "qgsrastertransparency.h"
#include "qgsrendercontext.h"
#include "qgsrenderedlayerstatistics.h"
#include "qgsruntimeprofiler.h"
#include "qgsthreadingutils.h"
#include "qgsunittypes.h"
 
#include <QElapsedTimer>
#include <QPointer>
#include <QString>
#include <QThread>
 
#include "moc_qgsrasterlayerrenderer.cpp"
 
using namespace Qt::StringLiterals;

 
QgsRasterLayerRendererFeedback::QgsRasterLayerRendererFeedback( QgsRasterLayerRenderer *r )
  : mR( r )
{
  setRenderPartialOutput( r->renderContext()->testFlag( Qgis::RenderContextFlag::RenderPartialOutput ) );
}
 
void QgsRasterLayerRendererFeedback::onNewData()
{
  if ( !renderPartialOutput() )
    return; // we were not asked for partial renders and we may not have a temporary image for overwriting...
 
  // update only once upon a time
  // (preview itself takes some time)
  if ( mLastPreview.isValid() && mLastPreview.msecsTo( QTime::currentTime() ) < mMinimalPreviewInterval )
    return;
 
  // TODO: update only the area that got new data
 
  QgsDebugMsgLevel( u"new raster preview! %1"_s.arg( mLastPreview.msecsTo( QTime::currentTime() ) ), 3 );
  QElapsedTimer t;
  t.start();
  QgsRasterBlockFeedback feedback;
  feedback.setPreviewOnly( true );
  feedback.setRenderPartialOutput( true );
  QgsRasterIterator iterator( mR->mPipe->last() );
  QgsRasterDrawer drawer( &iterator );
  drawer.draw( *( mR->renderContext() ), mR->mRasterViewPort, &feedback );
  mR->mReadyToCompose = true;
  QgsDebugMsgLevel( u"total raster preview time: %1 ms"_s.arg( t.elapsed() ), 3 );
  mLastPreview = QTime::currentTime();
}

QgsRasterLayerRenderer::QgsRasterLayerRenderer( QgsRasterLayer *layer, QgsRenderContext &rendererContext )
  : QgsMapLayerRenderer( layer->id(), &rendererContext )
  , mLayerName( layer->name() )
  , mLayerOpacity( layer->opacity() )
  , mProviderCapabilities( layer->dataProvider()->providerCapabilities() )
  , mInterfaceCapabilities( layer->dataProvider()->capabilities() )
  , mFeedback( new QgsRasterLayerRendererFeedback( this ) )
  , mEnableProfile( rendererContext.flags() & Qgis::RenderContextFlag::RecordProfile )
{
  mReadyToCompose = false;
 
  QElapsedTimer timer;
  timer.start();
 
  QgsMapToPixel mapToPixel = rendererContext.mapToPixel();
  if ( rendererContext.mapToPixel().mapRotation() )
  {
    // unset rotation for the sake of local computations.
    // Rotation will be handled by QPainter later
    // TODO: provide a method of QgsMapToPixel to fetch map center
    //       in geographical units
    const QgsPointXY center = mapToPixel.toMapCoordinates( static_cast<int>( mapToPixel.mapWidth() / 2.0 ), static_cast<int>( mapToPixel.mapHeight() / 2.0 ) );
    mapToPixel.setMapRotation( 0, center.x(), center.y() );
  }
 
  QgsRectangle viewExtentInMapCrs;
  QgsRectangle layerExtentInMapCrs;
 
  if ( rendererContext.coordinateTransform().isValid() )
  {
    QgsDebugMsgLevel( u"coordinateTransform set -> project extents."_s, 4 );
    if ( rendererContext.extent().xMinimum() == std::numeric_limits<double>::lowest()
         && rendererContext.extent().yMinimum() == std::numeric_limits<double>::lowest()
         && rendererContext.extent().xMaximum() == std::numeric_limits<double>::max()
         && rendererContext.extent().yMaximum() == std::numeric_limits<double>::max() )
    {
      // We get in this situation if the view CRS is geographical and the
      // extent goes beyond -180,-90,180,90. To avoid reprojection issues to the
      // layer CRS, then this dummy extent is returned by QgsMapRendererJob::reprojectToLayerExtent()
      // Don't try to reproject it now to view extent as this would return
      // a null rectangle.
      viewExtentInMapCrs = rendererContext.extent();
    }
    else
    {
      try
      {
        QgsCoordinateTransform ct = rendererContext.coordinateTransform();
        ct.setBallparkTransformsAreAppropriate( true );
        viewExtentInMapCrs = rendererContext.mapExtent();
      }
      catch ( QgsCsException &cs )
      {
        QgsMessageLog::logMessage( QObject::tr( "Could not reproject view extent: %1" ).arg( cs.what() ), QObject::tr( "Raster" ) );
        viewExtentInMapCrs.setNull();
      }
    }
 
    try
    {
      QgsCoordinateTransform ct = rendererContext.coordinateTransform();
      ct.setBallparkTransformsAreAppropriate( true );
      layerExtentInMapCrs = ct.transformBoundingBox( layer->extent() );
    }
    catch ( QgsCsException &cs )
    {
      QgsMessageLog::logMessage( QObject::tr( "Could not reproject layer extent: %1" ).arg( cs.what() ), QObject::tr( "Raster" ) );
      layerExtentInMapCrs.setNull();
    }
  }
  else
  {
    QgsDebugMsgLevel( u"coordinateTransform not set"_s, 4 );
    viewExtentInMapCrs = rendererContext.extent();
    layerExtentInMapCrs = layer->extent();
  }
 
  // clip raster extent to view extent
  QgsRectangle visibleExtentOfRasterInMapCrs = layer->ignoreExtents() ? viewExtentInMapCrs : viewExtentInMapCrs.intersect( layerExtentInMapCrs );
  if ( visibleExtentOfRasterInMapCrs.isEmpty() )
  {
    if ( rendererContext.coordinateTransform().isShortCircuited() )
    {
      // no point doing the fallback logic, we aren't using transforms
 
      QgsDebugMsgLevel( u"draw request outside view extent."_s, 2 );
      // nothing to do
      return;
    }
 
    // fallback logic, try the inverse operation
    // eg if we are trying to view a global layer in a small area local projection, then the transform of the
    // map extent to the layer's local projection probably failed.
    // So we'll try the opposite logic here, and see if we can first calculate the visible region in the
    // layers extent, and then transform that small region back to the map's crs
    try
    {
      QgsCoordinateTransform layerToMapTransform = rendererContext.coordinateTransform();
      layerToMapTransform.setBallparkTransformsAreAppropriate( true );
      const QgsRectangle viewExtentInLayerCrs = layerToMapTransform.transformBoundingBox( rendererContext.mapExtent(), Qgis::TransformDirection::Reverse );
 
      const QgsRectangle visibleExtentInLayerCrs = layer->ignoreExtents() ? viewExtentInLayerCrs : viewExtentInLayerCrs.intersect( layer->extent() );
      if ( visibleExtentInLayerCrs.isEmpty() )
      {
        QgsDebugMsgLevel( u"draw request outside view extent."_s, 2 );
        // nothing to do
        return;
      }
 
      visibleExtentOfRasterInMapCrs = layerToMapTransform.transformBoundingBox( visibleExtentInLayerCrs );
    }
    catch ( QgsCsException &cs )
    {
      QgsMessageLog::logMessage( QObject::tr( "Could not reproject layer extent: %1" ).arg( cs.what() ), QObject::tr( "Raster" ) );
      return;
    }
  }
 
  if ( visibleExtentOfRasterInMapCrs.isEmpty() )
  {
    QgsDebugMsgLevel( u"draw request outside view extent."_s, 2 );
    // nothing to do
    return;
  }
 
  QgsDebugMsgLevel( "map extent in layer crs is " + rendererContext.extent().toString(), 4 );
  QgsDebugMsgLevel( "map extent in map crs is " + viewExtentInMapCrs.toString(), 4 );
  QgsDebugMsgLevel( "layer extent in map crs is " + layerExtentInMapCrs.toString(), 4 );
  QgsDebugMsgLevel( "visible extent of raster in map crs is " + visibleExtentOfRasterInMapCrs.toString(), 4 );
 
  //
  // The first thing we do is set up the QgsRasterViewPort. This struct stores all the settings
  // relating to the size (in pixels and coordinate system units) of the raster part that is
  // in view in the map window. It also stores the origin.
  //
  //this is not a class level member because every time the user pans or zooms
  //the contents of the rasterViewPort will change
  mRasterViewPort = new QgsRasterViewPort();
 
  mRasterViewPort->mDrawnExtent = visibleExtentOfRasterInMapCrs;
  if ( rendererContext.coordinateTransform().isValid() )
  {
    mRasterViewPort->mSrcCRS = layer->crs();
    mRasterViewPort->mDestCRS = rendererContext.coordinateTransform().destinationCrs();
    mRasterViewPort->mTransformContext = rendererContext.transformContext();
  }
  else
  {
    mRasterViewPort->mSrcCRS = QgsCoordinateReferenceSystem();  // will be invalid
    mRasterViewPort->mDestCRS = QgsCoordinateReferenceSystem(); // will be invalid
  }
 
  // get dimensions of clipped raster image in device coordinate space (this is the size of the viewport)
  mRasterViewPort->mTopLeftPoint = mapToPixel.transform( visibleExtentOfRasterInMapCrs.xMinimum(), visibleExtentOfRasterInMapCrs.yMaximum() );
  mRasterViewPort->mBottomRightPoint = mapToPixel.transform( visibleExtentOfRasterInMapCrs.xMaximum(), visibleExtentOfRasterInMapCrs.yMinimum() );
 
  // align to output device grid, i.e. std::floor/ceil to integers
  // TODO: this should only be done if paint device is raster - screen, image
  // for other devices (pdf) it can have floating point origin
  // we could use floating point for raster devices as well, but respecting the
  // output device grid should make it more effective as the resampling is done in
  // the provider anyway
  mRasterViewPort->mTopLeftPoint.setX( std::floor( mRasterViewPort->mTopLeftPoint.x() ) );
  mRasterViewPort->mTopLeftPoint.setY( std::floor( mRasterViewPort->mTopLeftPoint.y() ) );
  mRasterViewPort->mBottomRightPoint.setX( std::ceil( mRasterViewPort->mBottomRightPoint.x() ) );
  mRasterViewPort->mBottomRightPoint.setY( std::ceil( mRasterViewPort->mBottomRightPoint.y() ) );
  // recalc myRasterExtent to aligned values
  visibleExtentOfRasterInMapCrs.set(
    mapToPixel.toMapCoordinates( mRasterViewPort->mTopLeftPoint.x(), mRasterViewPort->mBottomRightPoint.y() ),
    mapToPixel.toMapCoordinates( mRasterViewPort->mBottomRightPoint.x(), mRasterViewPort->mTopLeftPoint.y() )
  );
 
  //raster viewport top left / bottom right are already rounded to int
  mRasterViewPort->mWidth = static_cast<qgssize>( std::abs( mRasterViewPort->mBottomRightPoint.x() - mRasterViewPort->mTopLeftPoint.x() ) );
  mRasterViewPort->mHeight = static_cast<qgssize>( std::abs( mRasterViewPort->mBottomRightPoint.y() - mRasterViewPort->mTopLeftPoint.y() ) );
 
  double dpi = 25.4 * rendererContext.scaleFactor();
  if ( mProviderCapabilities & Qgis::RasterProviderCapability::DpiDependentData && rendererContext.dpiTarget() >= 0.0 )
  {
    const double dpiScaleFactor = rendererContext.dpiTarget() / dpi;
    mRasterViewPort->mWidth *= dpiScaleFactor;
    mRasterViewPort->mHeight *= dpiScaleFactor;
    dpi = rendererContext.dpiTarget();
  }
  else
  {
    rendererContext.setDpiTarget( -1.0 );
  }
 
  //the drawable area can start to get very very large when you get down displaying 2x2 or smaller, this is because
  //mapToPixel.mapUnitsPerPixel() is less then 1,
  //so we will just get the pixel data and then render these special cases differently in paintImageToCanvas()
 
  QgsDebugMsgLevel( u"mapUnitsPerPixel = %1"_s.arg( mapToPixel.mapUnitsPerPixel() ), 3 );
  QgsDebugMsgLevel( u"mWidth = %1"_s.arg( layer->width() ), 3 );
  QgsDebugMsgLevel( u"mHeight = %1"_s.arg( layer->height() ), 3 );
  QgsDebugMsgLevel( u"visibleExtentOfRasterInMapCrs.xMinimum() = %1"_s.arg( visibleExtentOfRasterInMapCrs.xMinimum() ), 3 );
  QgsDebugMsgLevel( u"visibleExtentOfRasterInMapCrs.xMaximum() = %1"_s.arg( visibleExtentOfRasterInMapCrs.xMaximum() ), 3 );
  QgsDebugMsgLevel( u"visibleExtentOfRasterInMapCrs.yMinimum() = %1"_s.arg( visibleExtentOfRasterInMapCrs.yMinimum() ), 3 );
  QgsDebugMsgLevel( u"visibleExtentOfRasterInMapCrs.yMaximum() = %1"_s.arg( visibleExtentOfRasterInMapCrs.yMaximum() ), 3 );
 
  QgsDebugMsgLevel( u"mTopLeftPoint.x() = %1"_s.arg( mRasterViewPort->mTopLeftPoint.x() ), 3 );
  QgsDebugMsgLevel( u"mBottomRightPoint.x() = %1"_s.arg( mRasterViewPort->mBottomRightPoint.x() ), 3 );
  QgsDebugMsgLevel( u"mTopLeftPoint.y() = %1"_s.arg( mRasterViewPort->mTopLeftPoint.y() ), 3 );
  QgsDebugMsgLevel( u"mBottomRightPoint.y() = %1"_s.arg( mRasterViewPort->mBottomRightPoint.y() ), 3 );
 
  QgsDebugMsgLevel( u"mWidth = %1"_s.arg( mRasterViewPort->mWidth ), 3 );
  QgsDebugMsgLevel( u"mHeight = %1"_s.arg( mRasterViewPort->mHeight ), 3 );
 
  // /\/\/\ - added to handle zoomed-in rasters
 
  // TODO R->mLastViewPort = *mRasterViewPort;
 
  // TODO: is it necessary? Probably WMS only?
  layer->dataProvider()->setDpi( std::floor( dpi * rendererContext.devicePixelRatio() ) );
 
  // copy the whole raster pipe!
  mPipe = std::make_unique<QgsRasterPipe>( *layer->pipe() );
 
  QObject::connect( mPipe->provider(), &QgsRasterDataProvider::statusChanged, layer, &QgsRasterLayer::statusChanged );
  QgsRasterRenderer *rasterRenderer = mPipe->renderer();
  if ( rasterRenderer && !( rendererContext.flags() & Qgis::RenderContextFlag::RenderPreviewJob ) && !( rendererContext.flags() & Qgis::RenderContextFlag::Render3DMap ) )
  {
    if ( rasterRenderer->needsRefresh( rendererContext.extent() ) )
    {
      QList<double> minValues;
      QList<double> maxValues;
      const QgsRasterMinMaxOrigin &minMaxOrigin = rasterRenderer->minMaxOrigin();
      for ( const int bandIdx : rasterRenderer->usesBands() )
      {
        double min;
        double max;
        layer->computeMinMax( bandIdx, minMaxOrigin, minMaxOrigin.limits(), rendererContext.extent(), static_cast<int>( QgsRasterLayer::SAMPLE_SIZE ), min, max );
        minValues.append( min );
        maxValues.append( max );
      }
 
      rasterRenderer->refresh( rendererContext.extent(), minValues, maxValues );
      QgsRenderedLayerStatistics *layerStatistics = new QgsRenderedLayerStatistics( layer->id(), minValues, maxValues );
      layerStatistics->setBoundingBox( rendererContext.extent() );
      appendRenderedItemDetails( layerStatistics );
    }
  }
 
  mPipe->evaluateDataDefinedProperties( rendererContext.expressionContext() );
 
  const QgsRasterLayerTemporalProperties *temporalProperties = qobject_cast< const QgsRasterLayerTemporalProperties * >( layer->temporalProperties() );
  const QgsRasterLayerElevationProperties *elevationProperties = qobject_cast<QgsRasterLayerElevationProperties *>( layer->elevationProperties() );
 
  if ( ( temporalProperties->isActive() && renderContext()->isTemporal() ) || ( elevationProperties->hasElevation() && !renderContext()->zRange().isInfinite() ) )
  {
    // temporal and/or elevation band filtering may be applicable
    bool matched = false;
    const int matchedBand = QgsRasterLayerUtils::renderedBandForElevationAndTemporalRange( layer, rendererContext.temporalRange(), rendererContext.zRange(), matched );
    if ( matched && matchedBand > 0 )
    {
      mPipe->renderer()->setInputBand( matchedBand );
    }
  }
 
  if ( temporalProperties->isActive() && renderContext()->isTemporal() )
  {
    switch ( temporalProperties->mode() )
    {
      case Qgis::RasterTemporalMode::FixedDateTime:
      case Qgis::RasterTemporalMode::FixedTemporalRange:
      case Qgis::RasterTemporalMode::RedrawLayerOnly:
      case Qgis::RasterTemporalMode::FixedRangePerBand:
        break;
 
      case Qgis::RasterTemporalMode::RepresentsTemporalValues:
        if ( mPipe->renderer()->usesBands().contains( temporalProperties->bandNumber() ) )
        {
          // if layer has elevation settings and we are only rendering a temporal range => we need to filter pixels by temporal values
          std::unique_ptr< QgsRasterTransparency > transparency;
          if ( const QgsRasterTransparency *rendererTransparency = mPipe->renderer()->rasterTransparency() )
            transparency = std::make_unique< QgsRasterTransparency >( *rendererTransparency );
          else
            transparency = std::make_unique< QgsRasterTransparency >();
 
          QVector<QgsRasterTransparency::TransparentSingleValuePixel> transparentPixels = transparency->transparentSingleValuePixelList();
 
          const QDateTime &offset = temporalProperties->temporalRepresentationOffset();
          const QgsInterval &scale = temporalProperties->temporalRepresentationScale();
          const double adjustedLower = static_cast< double >( offset.msecsTo( rendererContext.temporalRange().begin() ) )
                                       * QgsUnitTypes::fromUnitToUnitFactor( Qgis::TemporalUnit::Milliseconds, scale.originalUnit() )
                                       / scale.originalDuration();
          const double adjustedUpper = static_cast< double >( offset.msecsTo( rendererContext.temporalRange().end() ) )
                                       * QgsUnitTypes::fromUnitToUnitFactor( Qgis::TemporalUnit::Milliseconds, scale.originalUnit() )
                                       / scale.originalDuration();
          if ( !temporalProperties->accumulatePixels() )
          {
            transparentPixels.append( QgsRasterTransparency::TransparentSingleValuePixel( std::numeric_limits<double>::lowest(), adjustedLower, 0, true, !rendererContext.zRange().includeLower() ) );
          }
          transparentPixels.append( QgsRasterTransparency::TransparentSingleValuePixel( adjustedUpper, std::numeric_limits<double>::max(), 0, !rendererContext.zRange().includeUpper(), true ) );
 
          transparency->setTransparentSingleValuePixelList( transparentPixels );
          mPipe->renderer()->setRasterTransparency( transparency.release() );
        }
        break;
 
      case Qgis::RasterTemporalMode::TemporalRangeFromDataProvider:
        // in this mode we need to pass on the desired render temporal range to the data provider
        if ( QgsRasterDataProviderTemporalCapabilities *temporalCapabilities = mPipe->provider()->temporalCapabilities() )
        {
          temporalCapabilities->setRequestedTemporalRange( rendererContext.temporalRange() );
          temporalCapabilities->setIntervalHandlingMethod( temporalProperties->intervalHandlingMethod() );
        }
        break;
    }
  }
  else if ( QgsRasterDataProviderTemporalCapabilities *temporalCapabilities = mPipe->provider()->temporalCapabilities() )
  {
    temporalCapabilities->setRequestedTemporalRange( QgsDateTimeRange() );
    temporalCapabilities->setIntervalHandlingMethod( temporalProperties->intervalHandlingMethod() );
  }
 
  mClippingRegions = QgsMapClippingUtils::collectClippingRegionsForLayer( *renderContext(), layer );
 
  if ( elevationProperties && elevationProperties->hasElevation() )
  {
    mDrawElevationMap = true;
    mElevationScale = elevationProperties->zScale();
    mElevationOffset = elevationProperties->zOffset();
    mElevationBand = elevationProperties->bandNumber();
 
    if ( !rendererContext.zRange().isInfinite() )
    {
      // NOLINTBEGIN(bugprone-branch-clone)
      switch ( elevationProperties->mode() )
      {
        case Qgis::RasterElevationMode::FixedElevationRange:
          // don't need to handle anything here -- the layer renderer will never be created if the
          // render context range doesn't match the layer's fixed elevation range
          break;
 
        case Qgis::RasterElevationMode::FixedRangePerBand:
        case Qgis::RasterElevationMode::DynamicRangePerBand:
          // temporal/elevation band based filtering was already handled earlier in this method
          break;
 
        case Qgis::RasterElevationMode::RepresentsElevationSurface:
        {
          if ( mPipe->renderer()->usesBands().contains( mElevationBand ) )
          {
            // if layer has elevation settings and we are only rendering a slice of z values => we need to filter pixels by elevation
            if ( mPipe->renderer()->flags() & Qgis::RasterRendererFlag::UseNoDataForOutOfRangePixels )
            {
              std::unique_ptr< QgsRasterNuller> nuller;
              if ( const QgsRasterNuller *existingNuller = mPipe->nuller() )
                nuller.reset( existingNuller->clone() );
              else
                nuller = std::make_unique< QgsRasterNuller >();
 
              // account for z offset/zscale by reversing these calculations, so that we get the z range in
              // raw pixel values
              QgsRasterRangeList nullRanges;
              const double adjustedLower = ( rendererContext.zRange().lower() - mElevationOffset ) / mElevationScale;
              const double adjustedUpper = ( rendererContext.zRange().upper() - mElevationOffset ) / mElevationScale;
              nullRanges.append(
                QgsRasterRange( std::numeric_limits<double>::lowest(), adjustedLower, rendererContext.zRange().includeLower() ? QgsRasterRange::BoundsType::IncludeMin : QgsRasterRange::BoundsType::IncludeMinAndMax )
              );
              nullRanges.append(
                QgsRasterRange( adjustedUpper, std::numeric_limits<double>::max(), rendererContext.zRange().includeUpper() ? QgsRasterRange::BoundsType::IncludeMax : QgsRasterRange::BoundsType::IncludeMinAndMax )
              );
              nuller->setOutputNoDataValue( mElevationBand, static_cast< int >( adjustedLower - 1 ) );
              nuller->setNoData( mElevationBand, nullRanges );
 
              if ( !mPipe->insert( 1, nuller.release() ) )
              {
                QgsDebugError( u"Cannot set pipe nuller"_s );
              }
            }
            else
            {
              std::unique_ptr< QgsRasterTransparency > transparency;
              if ( const QgsRasterTransparency *rendererTransparency = mPipe->renderer()->rasterTransparency() )
                transparency = std::make_unique< QgsRasterTransparency >( *rendererTransparency );
              else
                transparency = std::make_unique< QgsRasterTransparency >();
 
              QVector<QgsRasterTransparency::TransparentSingleValuePixel> transparentPixels = transparency->transparentSingleValuePixelList();
 
              // account for z offset/zscale by reversing these calculations, so that we get the z range in
              // raw pixel values
              const double adjustedLower = ( rendererContext.zRange().lower() - mElevationOffset ) / mElevationScale;
              const double adjustedUpper = ( rendererContext.zRange().upper() - mElevationOffset ) / mElevationScale;
              transparentPixels.append( QgsRasterTransparency::TransparentSingleValuePixel( std::numeric_limits<double>::lowest(), adjustedLower, 0, true, !rendererContext.zRange().includeLower() ) );
              transparentPixels.append( QgsRasterTransparency::TransparentSingleValuePixel( adjustedUpper, std::numeric_limits<double>::max(), 0, !rendererContext.zRange().includeUpper(), true ) );
 
              transparency->setTransparentSingleValuePixelList( transparentPixels );
              mPipe->renderer()->setRasterTransparency( transparency.release() );
            }
          }
          break;
        }
      }
      // NOLINTEND(bugprone-branch-clone)
    }
  }
 
  prepareLabeling( layer );
 
  mFeedback->setRenderContext( rendererContext );
 
  mPipe->moveToThread( nullptr );
 
  mPreparationTime = timer.elapsed();
}
 
QgsRasterLayerRenderer::~QgsRasterLayerRenderer()
{
  delete mFeedback;
 
  delete mRasterViewPort;
}
 
bool QgsRasterLayerRenderer::render()
{
  QgsScopedThreadName threadName( u"render:%1"_s.arg( mLayerName ) );
 
  std::unique_ptr< QgsScopedRuntimeProfile > profile;
  if ( mEnableProfile )
  {
    profile = std::make_unique< QgsScopedRuntimeProfile >( mLayerName, u"rendering"_s, layerId() );
    if ( mPreparationTime > 0 )
      QgsApplication::profiler()->record( QObject::tr( "Create renderer" ), mPreparationTime / 1000.0, u"rendering"_s );
  }
 
  // Skip rendering of out of view tiles (xyz)
  if ( !mRasterViewPort || ( renderContext()->testFlag( Qgis::RenderContextFlag::RenderPreviewJob ) && !( mInterfaceCapabilities & Qgis::RasterInterfaceCapability::Prefetch ) ) )
    return true;
 
  mPipe->moveToThread( QThread::currentThread() );
 
  QElapsedTimer time;
  time.start();
 
  std::unique_ptr< QgsScopedRuntimeProfile > preparingProfile;
  if ( mEnableProfile )
  {
    preparingProfile = std::make_unique< QgsScopedRuntimeProfile >( QObject::tr( "Preparing render" ), u"rendering"_s );
  }
 
  //
  //
  // The goal here is to make as many decisions as possible early on (outside of the rendering loop)
  // so that we can maximise performance of the rendering process. So now we check which drawing
  // procedure to use :
  //
 
  const QgsScopedQPainterState painterSate( renderContext()->painter() );
 
  // important -- disable SmoothPixmapTransform and Antialiasing for raster layer renders. We want individual pixels to be clearly defined!
  renderContext()->painter()->setRenderHint( QPainter::SmoothPixmapTransform, false );
  bool antialiasEnabled = false;
  if ( !mClippingRegions.empty() )
  {
    bool needsPainterClipPath = false;
    const QPainterPath path = QgsMapClippingUtils::calculatePainterClipRegion( mClippingRegions, *renderContext(), Qgis::LayerType::Raster, needsPainterClipPath );
    if ( needsPainterClipPath )
    {
      renderContext()->painter()->setClipPath( path, Qt::IntersectClip );
      // this is an exception to the earlier flag disabling antialiasing!
      // If we are rendering with clip paths, we WANT the boundaries of the raster to be
      // rendered using antialiasing along the clip path, or things are very ugly! And since we can't
      // selectively antialias JUST the clip path boundary, we fallback to antialiasing the
      // whole raster render
      antialiasEnabled = true;
    }
  }
  renderContext()->painter()->setRenderHint( QPainter::Antialiasing, antialiasEnabled );
 
  QgsRasterProjector *projector = mPipe->projector();
  bool restoreOldResamplingStage = false;
  const Qgis::RasterResamplingStage oldResamplingState = mPipe->resamplingStage();
 
  // TODO add a method to interface to get provider and get provider
  // params in QgsRasterProjector
  if ( projector )
  {
    // Force provider resampling if reprojection is needed
    if ( ( mPipe->provider()->providerCapabilities() & Qgis::RasterProviderCapability::ProviderHintCanPerformProviderResampling )
         && mRasterViewPort->mSrcCRS != mRasterViewPort->mDestCRS
         && oldResamplingState != Qgis::RasterResamplingStage::Provider )
    {
      restoreOldResamplingStage = true;
      mPipe->setResamplingStage( Qgis::RasterResamplingStage::Provider );
    }
    projector->setCrs( mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS, mRasterViewPort->mTransformContext );
  }
 
  preparingProfile.reset();
  std::unique_ptr< QgsScopedRuntimeProfile > renderingProfile;
  if ( mEnableProfile )
  {
    renderingProfile = std::make_unique< QgsScopedRuntimeProfile >( QObject::tr( "Rendering" ), u"rendering"_s );
  }
 
  // Drawer to pipe?
  QgsRasterIterator iterator( mPipe->last() );
 
  // Get the maximum tile size from the provider and set it as the maximum tile size for the iterator
  if ( QgsRasterDataProvider *provider = mPipe->provider() )
  {
    const QSize maxTileSize { provider->maximumTileSize() };
    iterator.setMaximumTileWidth( maxTileSize.width() );
    iterator.setMaximumTileHeight( maxTileSize.height() );
  }
 
  QgsRasterDrawer drawer( &iterator );
  drawer.draw( *( renderContext() ), mRasterViewPort, mFeedback );
 
  if ( mDrawElevationMap )
    drawElevationMap();
 
  renderingProfile.reset();
 
  if ( mLabelProvider && !renderContext()->renderingStopped() )
  {
    std::unique_ptr< QgsScopedRuntimeProfile > labelingProfile;
    if ( mEnableProfile )
    {
      labelingProfile = std::make_unique< QgsScopedRuntimeProfile >( QObject::tr( "Labeling" ), u"rendering"_s );
    }
    drawLabeling();
  }
 
  if ( restoreOldResamplingStage )
  {
    mPipe->setResamplingStage( oldResamplingState );
  }
 
  const QStringList errors = mFeedback->errors();
  for ( const QString &error : errors )
  {
    mErrors.append( error );
  }
 
  QgsDebugMsgLevel( u"total raster draw time (ms):     %1"_s.arg( time.elapsed(), 5 ), 4 );
  mReadyToCompose = true;
 
  mPipe->moveToThread( nullptr );
 
  return !mFeedback->isCanceled();
}
 
QgsFeedback *QgsRasterLayerRenderer::feedback() const
{
  return mFeedback;
}
 
bool QgsRasterLayerRenderer::forceRasterRender() const
{
  if ( !mRasterViewPort || !mPipe )
    return false; // this layer is not going to get rendered
 
  // preview of intermediate raster rendering results requires a temporary output image
  if ( renderContext()->testFlag( Qgis::RenderContextFlag::RenderPartialOutput ) )
    return true;
 
  if ( QgsRasterRenderer *renderer = mPipe->renderer() )
  {
    if ( !( renderer->flags() & Qgis::RasterRendererFlag::InternalLayerOpacityHandling ) )
    {
      switch ( renderContext()->rasterizedRenderingPolicy() )
      {
        case Qgis::RasterizedRenderingPolicy::Default:
        case Qgis::RasterizedRenderingPolicy::PreferVector:
          break;
 
        case Qgis::RasterizedRenderingPolicy::ForceVector:
          return false;
      }
 
      if ( !qgsDoubleNear( mLayerOpacity, 1.0 ) )
        return true;
    }
  }
 
  return false;
}
 
void QgsRasterLayerRenderer::prepareLabeling( QgsRasterLayer *layer )
{
  QgsRenderContext &context = *renderContext();
 
  if ( QgsLabelingEngine *engine2 = context.labelingEngine() )
  {
    if ( QgsAbstractRasterLayerLabeling *labeling = layer->labeling() )
    {
      if ( layer->labelsEnabled() && labeling->isInScaleRange( context.rendererScale() ) )
      {
        std::unique_ptr< QgsRasterLayerLabelProvider > provider = labeling->provider( layer );
        if ( provider )
        {
          // engine takes ownership
          mLabelProvider = provider.release();
          mLabelProvider->startRender( context );
          engine2->addProvider( mLabelProvider );
        }
      }
    }
  }
}
 
void QgsRasterLayerRenderer::drawLabeling()
{
  if ( mLabelProvider )
    mLabelProvider->generateLabels( *renderContext(), mPipe.get(), mRasterViewPort, mFeedback );
}
 
void QgsRasterLayerRenderer::drawElevationMap()
{
  QgsRasterDataProvider *dataProvider = mPipe->provider();
  if ( renderContext()->elevationMap() && dataProvider )
  {
    double dpiScalefactor;
 
    if ( renderContext()->dpiTarget() >= 0.0 )
      dpiScalefactor = renderContext()->dpiTarget() / ( renderContext()->scaleFactor() * 25.4 );
    else
      dpiScalefactor = 1.0;
 
    int outputWidth = static_cast<int>( static_cast<double>( mRasterViewPort->mWidth ) / dpiScalefactor * renderContext()->devicePixelRatio() );
    int outputHeight = static_cast<int>( static_cast<double>( mRasterViewPort->mHeight ) / dpiScalefactor * renderContext()->devicePixelRatio() );
 
    QSize viewSize = renderContext()->deviceOutputSize();
    int viewWidth = static_cast<int>( viewSize.width() / dpiScalefactor );
    int viewHeight = static_cast<int>( viewSize.height() / dpiScalefactor );
 
    bool canRenderElevation = false;
    std::unique_ptr<QgsRasterBlock> elevationBlock;
    if ( mRasterViewPort->mSrcCRS == mRasterViewPort->mDestCRS )
    {
      elevationBlock.reset( dataProvider->block( mElevationBand, mRasterViewPort->mDrawnExtent, outputWidth, outputHeight, mFeedback ) );
      canRenderElevation = true;
    }
    else
    {
      // Destinaton CRS is different from the source CRS.
      // Using the raster projector lead to have big artifacts when rendering the elevation map.
      // To get a smoother elevation map, we use GDAL resampling with coordinates transform
      QgsRectangle viewExtentInLayerCoordinate = renderContext()->extent();
 
      // If view extent is infinite, we use the data provider extent
      if ( viewExtentInLayerCoordinate.xMinimum() == std::numeric_limits<double>::lowest()
           && viewExtentInLayerCoordinate.yMinimum() == std::numeric_limits<double>::lowest()
           && viewExtentInLayerCoordinate.xMaximum() == std::numeric_limits<double>::max()
           && viewExtentInLayerCoordinate.yMaximum() == std::numeric_limits<double>::max() )
      {
        viewExtentInLayerCoordinate = dataProvider->extent();
      }
 
      double xLayerResol = viewExtentInLayerCoordinate.width() / static_cast<double>( viewWidth );
      double yLayerResol = viewExtentInLayerCoordinate.height() / static_cast<double>( viewHeight );
 
      double overSampling = 1;
      if ( mPipe->resampleFilter() )
        overSampling = mPipe->resampleFilter()->maxOversampling();
 
      if ( dataProvider->capabilities() & Qgis::RasterInterfaceCapability::Size )
      {
        // If the dataprovider has size capability, we calculate the requested resolution to provider
        double providerXResol = dataProvider->extent().width() / dataProvider->xSize();
        double providerYResol = dataProvider->extent().height() / dataProvider->ySize();
        overSampling = ( xLayerResol / providerXResol + yLayerResol / providerYResol ) / 2;
      }
 
      GDALResampleAlg alg;
      if ( overSampling > 1 )
        alg = QgsGdalUtils::gdalResamplingAlgorithm( dataProvider->zoomedOutResamplingMethod() );
      else
        alg = QgsGdalUtils::gdalResamplingAlgorithm( dataProvider->zoomedInResamplingMethod() );
 
      Qgis::DataType dataType = dataProvider->dataType( mElevationBand );
 
      if ( dataType != Qgis::DataType::UnknownDataType ) // resampling data by GDAL is not compatible with unknown data type
      {
        // we need extra pixels on border to avoid effect border with resampling (at least 2 pixels band for cubic alg)
        int sourceWidth = viewWidth + 4;
        int sourceHeight = viewHeight + 4;
        viewExtentInLayerCoordinate = QgsRectangle(
          viewExtentInLayerCoordinate.xMinimum() - xLayerResol * 2,
          viewExtentInLayerCoordinate.yMinimum() - yLayerResol * 2,
          viewExtentInLayerCoordinate.xMaximum() + xLayerResol * 2,
          viewExtentInLayerCoordinate.yMaximum() + yLayerResol * 2
        );
 
        // Now we can do the resampling
        std::unique_ptr<QgsRasterBlock> sourcedata( dataProvider->block( mElevationBand, viewExtentInLayerCoordinate, sourceWidth, sourceHeight, mFeedback ) );
        gdal::dataset_unique_ptr gdalDsInput = QgsGdalUtils::blockToSingleBandMemoryDataset( viewExtentInLayerCoordinate, sourcedata.get() );
 
 
        elevationBlock = std::make_unique<QgsRasterBlock>( dataType, outputWidth, outputHeight );
 
        elevationBlock->setNoDataValue( dataProvider->sourceNoDataValue( mElevationBand ) );
 
        gdal::dataset_unique_ptr gdalDsOutput = QgsGdalUtils::blockToSingleBandMemoryDataset( mRasterViewPort->mDrawnExtent, elevationBlock.get() );
 
        // For coordinate transformation, we try to obtain a coordinate operation string from the transform context.
        // Depending of the CRS, if we can't we use GDAL transformation directly from the source and destination CRS
        QString coordinateOperation;
        const QgsCoordinateTransformContext &transformContext = renderContext()->transformContext();
        if ( transformContext.mustReverseCoordinateOperation( mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS ) )
          coordinateOperation = transformContext.calculateCoordinateOperation( mRasterViewPort->mDestCRS, mRasterViewPort->mSrcCRS );
        else
          coordinateOperation = transformContext.calculateCoordinateOperation( mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS );
 
        if ( coordinateOperation.isEmpty() )
          canRenderElevation = QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), alg, mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS );
        else
          canRenderElevation = QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), alg, coordinateOperation.toUtf8().constData() );
      }
    }
 
    if ( canRenderElevation )
    {
      QPoint topLeft;
      if ( renderContext()->mapToPixel().mapRotation() != 0 )
      {
        // Now rendering elevation on the elevation map, we need to take care of rotation:
        // again a resampling but this time with a geotransform.
        const QgsMapToPixel &mtp = renderContext()->mapToPixel();
        QgsElevationMap *elevMap = renderContext()->elevationMap();
 
        int elevMapWidth = elevMap->rawElevationImage().width();
        int elevMapHeight = elevMap->rawElevationImage().height();
 
        int bottom = 0;
        int top = elevMapHeight;
        int left = elevMapWidth;
        int right = 0;
 
        QList<QgsPointXY> corners;
        corners
          << QgsPointXY( mRasterViewPort->mDrawnExtent.xMinimum(), mRasterViewPort->mDrawnExtent.yMinimum() )
          << QgsPointXY( mRasterViewPort->mDrawnExtent.xMaximum(), mRasterViewPort->mDrawnExtent.yMaximum() )
          << QgsPointXY( mRasterViewPort->mDrawnExtent.xMinimum(), mRasterViewPort->mDrawnExtent.yMaximum() )
          << QgsPointXY( mRasterViewPort->mDrawnExtent.xMaximum(), mRasterViewPort->mDrawnExtent.yMinimum() );
 
        for ( const QgsPointXY &corner : std::as_const( corners ) )
        {
          const QgsPointXY dpt = mtp.transform( corner );
          int x = static_cast<int>( std::round( dpt.x() ) );
          int y = static_cast<int>( std::round( dpt.y() ) );
 
          if ( x < left )
            left = x;
          if ( x > right )
            right = x;
          if ( y < top )
            top = y;
          if ( y > bottom )
            bottom = y;
        }
 
        const QgsPointXY origin = mtp.toMapCoordinates( left, top );
        double gridXSize = mtp.toMapCoordinates( right, top ).distance( origin );
        double gridYSize = mtp.toMapCoordinates( left, bottom ).distance( origin );
        double angleRad = renderContext()->mapToPixel().mapRotation() / 180 * M_PI;
 
        gdal::dataset_unique_ptr gdalDsInput = QgsGdalUtils::blockToSingleBandMemoryDataset( mRasterViewPort->mDrawnExtent, elevationBlock.get() );
 
        std::unique_ptr<QgsRasterBlock> rotatedElevationBlock
          = std::make_unique<QgsRasterBlock>( elevationBlock->dataType(), ( right - left ) * renderContext()->devicePixelRatio() + 1, ( bottom - top ) * renderContext()->devicePixelRatio() + 1 );
 
        rotatedElevationBlock->setNoDataValue( elevationBlock->noDataValue() );
 
        gdal::dataset_unique_ptr gdalDsOutput = QgsGdalUtils::blockToSingleBandMemoryDataset( angleRad, origin, gridXSize, gridYSize, rotatedElevationBlock.get() );
 
        if ( QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), QgsGdalUtils::gdalResamplingAlgorithm( dataProvider->zoomedInResamplingMethod() ), nullptr ) )
        {
          elevationBlock = std::move( rotatedElevationBlock );
        }
 
        topLeft = QPoint( left, top );
      }
      else
      {
        topLeft = mRasterViewPort->mTopLeftPoint.toQPointF().toPoint();
      }
 
      renderContext()->elevationMap()->fillWithRasterBlock( elevationBlock.get(), topLeft.y() * renderContext()->devicePixelRatio(), topLeft.x() * renderContext()->devicePixelRatio(), mElevationScale, mElevationOffset );
    }
  }
}