qgsrasterlayer.cpp

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

  mUserDefinedGrayMinMaxFlag = false;

  mRasterShader = new QgsRasterShader();

  if ( loadDefaultStyleFlag )
  {
    bool defaultLoadedFlag = false;
    loadDefaultStyle( defaultLoadedFlag );
    if ( defaultLoadedFlag )
    {
      return;
    }
  }

  // Initialise the affine transform matrix
  mGeoTransform[0] =  0;
  mGeoTransform[1] =  1;
  mGeoTransform[2] =  0;
  mGeoTransform[3] =  0;
  mGeoTransform[4] =  0;
  mGeoTransform[5] = -1;

  // set the layer name (uppercase first character)

  if ( ! baseName.isEmpty() )   // XXX shouldn't this happen in parent?
  {
    setLayerName(baseName);
  }

  // load the file if one specified
  if ( ! path.isEmpty() )
  {
    readFile( path ); // XXX check for failure?
  }

} // QgsRasterLayer ctor


QgsRasterLayer::~QgsRasterLayer()
{

  if (mProviderKey.isEmpty())
  {
    if ( mGdalBaseDataset )
    {
      GDALDereferenceDataset( mGdalBaseDataset );
      GDALClose(mGdalDataset);
    }
  }  
}



bool QgsRasterLayer::readFile( QString const & fileName )
{
  registerGdalDrivers();

  mGdalDataset = NULL;

  //open the dataset making sure we handle char encoding of locale properly
  mGdalBaseDataset = GDALOpen(QFile::encodeName(fileName).constData(), GA_ReadOnly);

  if (mGdalBaseDataset == NULL)
  {
    mValid = FALSE;
    return false;
  }

  // Store timestamp
  mLastModified = lastModified ( fileName );

  // Check if we need a warped VRT for this file. 
  if( (GDALGetGeoTransform(mGdalBaseDataset,mGeoTransform) == CE_None
       && (mGeoTransform[1] < 0.0
           || mGeoTransform[2] != 0.0
           || mGeoTransform[4] != 0.0
           || mGeoTransform[5] > 0.0))
       || GDALGetGCPCount(mGdalBaseDataset) > 0 )
  {
      QgsLogger::warning("Creating Warped VRT.");

      mGdalDataset = 
          GDALAutoCreateWarpedVRT( mGdalBaseDataset, NULL, NULL, 
                                   GRA_NearestNeighbour, 0.2, NULL );
      if( mGdalDataset == NULL )
      {
          QgsLogger::warning("Warped VRT Creation failed.");
          mGdalDataset = mGdalBaseDataset;
          GDALReferenceDataset( mGdalDataset );
      }
  }
  else
  {
      mGdalDataset = mGdalBaseDataset;
      GDALReferenceDataset( mGdalDataset );
  }

  //check f this file has pyramids
  GDALRasterBandH myGDALBand = GDALGetRasterBand( mGdalDataset, 1 ); //just use the first band
  if( GDALGetOverviewCount(myGDALBand) > 0 )
  {
    hasPyramidsFlag=true;
  }
  else
  {
    hasPyramidsFlag=false;
  }

  //populate the list of what pyramids exist
  buildRasterPyramidList();

  // Get the layer's projection info and set up the
  // QgsCoordinateTransform for this layer
  // NOTE: we must do this before getMetadata is called

  QgsDebugMsg("Raster initial SRS");
  mSRS->debugPrint();

  QString mySourceWKT = getProjectionWKT();

  QgsDebugMsg("--------------------------------------------------------------------------------------");
  QgsDebugMsg("QgsRasterLayer::readFile --- using wkt\n" + mySourceWKT);
  QgsDebugMsg("--------------------------------------------------------------------------------------");

  mSRS->createFromWkt(mySourceWKT);
  //get the project projection, defaulting to this layer's projection
  //if none exists....
  if (!mSRS->isValid())
  {
    mSRS->validate();
  }
  QgsDebugMsg("Raster determined to have the following SRS");
  mSRS->debugPrint();

  //set up the coordinat transform - in the case of raster this is mainly used to convert
  //the inverese projection of the map extents of the canvas when zzooming in etc. so
  //that they match the coordinate system of this layer
  QgsDebugMsg("Layer registry has " + QString::number(QgsMapLayerRegistry::instance()->count()) + "layers");

  getMetadata();

  // Use the affine transform to get geo coordinates for
  // the corners of the raster
  double myXMax = mGeoTransform[0] +
    GDALGetRasterXSize(mGdalDataset) * mGeoTransform[1] +
    GDALGetRasterYSize(mGdalDataset) * mGeoTransform[2];
  double myYMin = mGeoTransform[3] +
    GDALGetRasterXSize(mGdalDataset) * mGeoTransform[4] +
    GDALGetRasterYSize(mGdalDataset) * mGeoTransform[5];

  mLayerExtent.setXmax(myXMax);
  // The affine transform reduces to these values at the
  // top-left corner of the raster
  mLayerExtent.setXmin(mGeoTransform[0]);
  mLayerExtent.setYmax(mGeoTransform[3]);
  mLayerExtent.setYmin(myYMin);

  //
  // Set up the x and y dimensions of this raster layer
  //
  mRasterXDim = GDALGetRasterXSize(mGdalDataset);
  mRasterYDim = GDALGetRasterYSize(mGdalDataset);

  //
  // Determin the nodatavalue
  //
  mNoDataValue = -9999; //Standard default?
  mValidNoDataValue = false;
  int isValid = false;
  double myNoDataValue = GDALGetRasterNoDataValue(GDALGetRasterBand(mGdalDataset,1),&isValid);
  if(isValid)
  {
    mNoDataValue = myNoDataValue;
    mValidNoDataValue = true;
  }
  
  if(mValidNoDataValue)
  {
    mRasterTransparency.initializeTransparentPixelList(mNoDataValue, mNoDataValue, mNoDataValue);
    mRasterTransparency.initializeTransparentPixelList(mNoDataValue);
  }

  //initialise the raster band stats and contrast enhancement vector
  for (int i = 1; i <= GDALGetRasterCount(mGdalDataset); i++)
  {
    GDALRasterBandH myGdalBand = GDALGetRasterBand(mGdalDataset,i);
    QString myColorQString = GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(myGdalBand));
    QgsRasterBandStats myRasterBandStats;
    //myRasterBandStats.bandName = myColorQString ;
    myRasterBandStats.bandName=QString::number(i) + " : " + myColorQString;
    myRasterBandStats.bandNo = i;
    myRasterBandStats.statsGatheredFlag = false;
    myRasterBandStats.histogramVector = new QgsRasterBandStats::HistogramVector();
    // Read color table
    readColorTable ( myGdalBand, &(myRasterBandStats.colorTable) );

    mRasterStatsList.push_back(myRasterBandStats);
    
    //Build a new contrast enhancement for the band and store in list
    QgsContrastEnhancement myContrastEnhancement((QgsContrastEnhancement::QgsRasterDataType)GDALGetRasterDataType(myGdalBand));
    mContrastEnhancementList.append(myContrastEnhancement);
  }
  
  //defaults - Needs to be set after the Contrast list has been build
  setContrastEnhancementAlgorithm(QgsContrastEnhancement::STRETCH_TO_MINMAX);
  
  //decide what type of layer this is...
  //note that multiband images can have one or more 'undefindd' bands,
  //so we must do this check first!
  if ((GDALGetRasterCount(mGdalDataset) > 1))
  {
    rasterLayerType = MULTIBAND;
  }
  else if (hasBand("Palette")) //dont tr() this its a gdal word!
  {
    rasterLayerType = PALETTE;
  }
  else
  {
    rasterLayerType = GRAY_OR_UNDEFINED;
  }

  if (rasterLayerType == PALETTE)
  {
    mRedBandName = "Red"; // sensible default
    mGreenBandName = "Green"; // sensible default
    mBlueBandName = "Blue";// sensible default
    mTransparencyBandName = tr(QSTRING_NOT_SET); // sensible default
    mGrayBandName = tr(QSTRING_NOT_SET);  //sensible default
    drawingStyle = PALETTED_MULTI_BAND_COLOR; //sensible default
    setContrastEnhancementAlgorithm(QgsContrastEnhancement::NO_STRETCH);
  }
  else if (rasterLayerType == MULTIBAND)
  {
    //we know we have at least 2 layers...
    mRedBandName = getRasterBandName(1);  // sensible default
    mGreenBandName = getRasterBandName(2);  // sensible default
    //for the third layer we cant be sure so..
    if (GDALGetRasterCount(mGdalDataset) > 2)
    {
      mBlueBandName = getRasterBandName(3); // sensible default
    }
    else
    {
      mBlueBandName = tr(QSTRING_NOT_SET);  // sensible default
    }

    //Beacuse the transparent band can be set from a different layer defaulting to the fourth
    //is not a sensible default P.Ersts 2007-05-14
    /*
    if (mGdalDataset->GetRasterCount() > 3)
      mTransparencyBandName = getRasterBandName(4);
    else
    */
    mTransparencyBandName = tr(QSTRING_NOT_SET);

    mGrayBandName = tr(QSTRING_NOT_SET);  //sensible default
    drawingStyle = MULTI_BAND_COLOR;  //sensible default
  }
  else                        //GRAY_OR_UNDEFINED
  {
    //Disabled automatically generating stats to improve initial load speed.
    //getRasterBandStats(1);
    mRedBandName = tr(QSTRING_NOT_SET); //sensible default
    mGreenBandName = tr(QSTRING_NOT_SET); //sensible default
    mBlueBandName = tr(QSTRING_NOT_SET);  //sensible default
    mTransparencyBandName = tr(QSTRING_NOT_SET);  //sensible default
    drawingStyle = SINGLE_BAND_GRAY;  //sensible default
    mGrayBandName = getRasterBandName(1); // usually gdal will return gray or undefined  
  }

  //mark the layer as valid
  mValid=TRUE;
  return true;

} // QgsRasterLayer::readFile

QString QgsRasterLayer::getProjectionWKT() 
{ 
  QString myWKTString;
  QgsSpatialRefSys mySRS;   
  myWKTString=QString (GDALGetProjectionRef(mGdalDataset));
  mySRS.createFromWkt(myWKTString);
  if (!mySRS.isValid())
  {
    //try to get the gcp srs from the raster layer if available
    myWKTString=QString(GDALGetGCPProjection(mGdalDataset));

// What is the purpose of this piece of code?
// Sideeffects from validate()?
//    mySRS.createFromWkt(myWKTString);
//    if (!mySRS.isValid())
//    {
//      // use force and make SRS valid!
//      mySRS.validate();
//    }

  }

  return myWKTString;
}

void QgsRasterLayer::closeDataset()
{
  if ( !mValid  ) return;
  mValid = FALSE;

  GDALDereferenceDataset(mGdalBaseDataset);
  mGdalBaseDataset = NULL;

  GDALClose(mGdalDataset);
  mGdalDataset = NULL;

  hasPyramidsFlag=false;
  mPyramidList.clear();

  mRasterStatsList.clear();
}

bool QgsRasterLayer::update()
{
  QgsDebugMsg("QgsRasterLayer::update");

  if ( mLastModified < QgsRasterLayer::lastModified ( source() ) )
  {
    QgsDebugMsg("Outdated -> reload");
    closeDataset();
    return readFile ( source() );
  }
  return true;
}

QDateTime QgsRasterLayer::lastModified ( QString const & name )
{
  QgsDebugMsg("QgsRasterLayer::lastModified: " + name);
  QDateTime t;

  QFileInfo fi ( name );

  // Is it file?
  if ( !fi.exists() ) return t;

  t = fi.lastModified();

  // Check also color table for GRASS
  if ( name.contains( "cellhd" ) > 0 )
  { // most probably GRASS
    QString dir = fi.dirPath();
    QString map = fi.fileName();
    fi.setFile ( dir + "/../colr/" + map );

    if ( fi.exists() )
    {
      if ( fi.lastModified() > t ) t = fi.lastModified();
    }
  }

  // Check GRASS group members (bands)
  if ( name.contains( "group" ) > 0 )
  { // probably GRASS group
    fi.setFile ( name + "/REF" );

    if ( fi.exists() )
    {  // most probably GRASS group
      QFile f ( name + "/REF" );
      if ( f.open ( QIODevice::ReadOnly ) )
      {
        QString dir = fi.dirPath() + "/../../../";

        char buf[101];
        while ( f.readLine(buf,100) != -1 )
        {
          QString ln = QString(buf);
          QStringList sl = QStringList::split ( ' ', ln.stripWhiteSpace() );
          QString map = sl.first();
          sl.pop_front();
          QString mapset = sl.first();

          // header
          fi.setFile ( dir + mapset + "/cellhd/" +  map );
          if ( fi.exists() )
          {
            if ( fi.lastModified() > t ) t = fi.lastModified();
          }

          // color
          fi.setFile ( dir + mapset + "/colr/" +  map );
          if ( fi.exists() )
          {
            if ( fi.lastModified() > t ) t = fi.lastModified();