GLIMSView Source Code Documentation

gdalfrmt.cpp

Go to the documentation of this file.

#include "gdalfrmt.h"
#include <iostream>
#include "TraceLog.h"

GDALFrmt::GDALFrmt( ) {
  GDALAllRegister();
  mDataSet = NULL;
}

GDALFrmt::~GDALFrmt( ) {
  if ( mDataSet )
    GDALClose( mDataSet );
}

void* GDALFrmt::getBand( int iband,
                         int x1,
                         int y1,
                         int w1,
                         int h1,
                         int w2,
                         int h2,
                         ScaleType ) const {
  int dSize      = 1;

  if ( mBandSet[iband].type == ImageFormat::UINT16 ||
       mBandSet[iband].type == ImageFormat::INT16 )
    dSize = 2;

  unsigned char *data = new unsigned char[ w2 * h2 * dSize ];

  GDALRasterBand *poBand = mDataSet->GetRasterBand( iband +1 );
  if ( !poBand ) {
    delete [] data;
    return NULL;
  }

  GDALDataType dt;
  switch ( mBandSet[iband].type ) {
    case ImageFormat::UINT8:
      dt = GDT_Byte;
      break;

    case ImageFormat::UINT16:
      dt = GDT_UInt16;
      break;

    case ImageFormat::INT16:
      dt = GDT_Int16;
      break;

    default:
      return data;
  }

  poBand->RasterIO( GF_Read, x1, y1, w1, h1, data, w2, h2, dt, 0, 0 );
  return(void*)data;
}

void GDALFrmt::getLL( int,
                      double &i,
                      double &j ) {
  double set[6];
  OGRSpatialReference poORIG, *poLL;
  const char *pszUTM_WKT;

  pszUTM_WKT = mDataSet->GetProjectionRef();

  if ( poORIG.importFromWkt( (char**)&pszUTM_WKT ) != OGRERR_NONE )
    return;

  if ( mDataSet->GetGeoTransform( set ) == CE_Failure )
    return;

  double dfGeorefX = set[0] + (set[1] * i) + (set[2] * j);
  double dfGeorefY = set[3] + (set[4] * i) + (set[5] * j);

  poLL = poORIG.CloneGeogCS();

  OGRCoordinateTransformation *poTransform = NULL;
  poTransform = OGRCreateCoordinateTransformation( &poORIG, poLL );

  if ( !poTransform )
    return;

  poTransform->Transform( 1, &dfGeorefX, &dfGeorefY );

  i = dfGeorefX;
  j = dfGeorefY;

  delete poLL;
  delete poTransform;
}

void GDALFrmt::getXY( int,
                      double &i,
                      double &j ) {
  OGRSpatialReference poORIG, *poLL;
  double set[6];
  const char *pszUTM_WKT;


  pszUTM_WKT = mDataSet->GetProjectionRef();

  if ( poORIG.importFromWkt( (char**)&pszUTM_WKT ) != OGRERR_NONE )
    return;

  if ( mDataSet->GetGeoTransform( set ) == CE_Failure )
    return;

  poLL = poORIG.CloneGeogCS();
  OGRCoordinateTransformation *poTransform = NULL;
  poTransform = OGRCreateCoordinateTransformation( poLL, &poORIG );
  if ( !poTransform ) {
    std::cerr << "GV_GDAL: Failed to get geotransform\n";
    return;
  }

  poTransform->Transform( 1, &i, &j );

  // set[0] + (set[1] * i) + (set[2] * j);
  // set[3] + (set[4] * i) + (set[5] * j);
  i = (i - set[0]) / set[1];
  j = (j - set[3]) / set[5];

  delete poLL;
  // delete poORIG;
  delete poTransform;
}

bool GDALFrmt::open( std::string fname ) {
  int     gotMin, gotMax;
  int     dims[2];
  int     nBand;
  double  minmax[2];
  char    buf[150];
  int     i=0;
  double  spattr[6];
  bool    isLL = false;
  const   char *pszUTM_WKT;

  OGRSpatialReference poOGR;
  GDALRasterBand *poBand;
  mFName = fname;

  mDataSet = (GDALDataset*)GDALOpen( fname.c_str(), GA_ReadOnly );
  if ( !mDataSet )
    return false;

  nBand = mDataSet->GetRasterCount();
  mDataSet->GetGeoTransform( spattr );
  pszUTM_WKT = mDataSet->GetProjectionRef();
  if ( poOGR.importFromWkt( (char**)&pszUTM_WKT ) == OGRERR_NONE )
    isLL = poOGR.IsGeographic();

  for ( i=0; i < nBand; i++ ) {
    poBand = mDataSet->GetRasterBand( i+1 );
    if ( !poBand )
      continue;

    dims[0] = poBand->GetXSize();
    dims[1] = poBand->GetYSize();

    minmax[0] = poBand->GetMinimum( &gotMin );
    minmax[1] = poBand->GetMaximum( &gotMax );
    if ( !(gotMin && gotMax) )
      GDALComputeRasterMinMax( poBand, true, minmax );

    sprintf( buf, "Band %d", i );

    ImageFormat::DataType dt;
    switch ( poBand->GetRasterDataType() ) {
      case GDT_Byte:
        dt = ImageFormat::UINT8;
        break;

      case GDT_UInt16:
        dt = ImageFormat::UINT16;
        break;

      case GDT_Int16:
        dt = ImageFormat::INT16;
        break;

      default:
        dt = ImageFormat::NOT_DEFINED;
        break;
    }

    BandInf b;
    b.width     = dims[0];
    b.height    = dims[1];
    b.type      = dt;
    b.name      = buf;
    b.min       = minmax[0];
    b.max       = minmax[1];
    b.ndims     = 2;
    b.res       = spattr[5];
    if ( b.res < 0 )
      b.res *= -1;

    if ( isLL )
      b.res *= 75000;
    mBandSet.push_back( b );
  }

  if ((int)mBandSet.size()) {
    double x = mBandSet[0].width / 2;
    double y = mBandSet[0].height / 2;
    getLL( 0, x, y );
    mImageInf.mCenterLat = y;
    mImageInf.mCenterLon = x;
  }
  mImageInf.mLocURL = fname ;

  return true;
}

---

Home | Search | Disclaimers & Privacy | Contact Us GLIMSView Maintainer: dsoltesz@usgs.gov