GLIMSView Source Code Documentation

hdfeosfrmt.cpp

Go to the documentation of this file.

#include "hdfeosfrmt.h"
#include <qmessagebox.h>
#include "TraceLog.h"

HdfEosFrmt::~HdfEosFrmt( ) {
  SWclose( mSWfid );
  for ( int i=0; i < (int)mGFList.size(); i++ ) {
    GeoField &gf = mGFList[i];
    switch ( gf.type ) {
      case DFNT_FLOAT32:
        delete [] (float*)gf.geomap;
        break;

      case DFNT_FLOAT64:
        delete [] (double*)gf.geomap;
        break;
    }
  }
}

bool HdfEosFrmt::open( std::string fname ) {
  char                    *list;
  std::vector<std::string> swlist;
  std::vector<std::string> fldlist;
  int32                    nret, nswath, bufsz, swid;
  int                      i;

  mFName = fname;
  if ( mFName == "" ) return false;

  mSWfid = SWopen( (char*)mFName.c_str(), DFACC_READ );
  if ( mSWfid == -1 ) return false;

  // GET NUMBER OF SWATHS AND LIST STRING SIZE
  nswath = SWinqswath( (char*)mFName.c_str(), NULL, &bufsz );
  if ( nswath < 1 || bufsz < 1 ) return false;

  list = new char[bufsz+2];
  list[bufsz+1] = (char)0;
  SWinqswath( (char*)mFName.c_str(), list, &bufsz );
  swlist = splitList( list, ',' );
  delete [] list;

  for ( i=0; i < nswath; i++ ) {
    int nfield;
    int32 *rank, *type, *dims;

    swid = SWattach( mSWfid, (char*)swlist[i].c_str() );
    if ( swid == -1 ) continue;  // UNABLE TO ATTACH

    parseGeoFlds( swid, swlist[i] );

    nret = SWnentries( swid, HDFE_NENTDFLD, &bufsz );
    if ( nret <= 0 ) continue;  // NO DATA FIELDS

    rank = new int32[nret];
    type = new int32[nret];
    list = new char[bufsz];
    nfield = SWinqdatafields( swid, list, rank, type );

    fldlist = splitList( list, ',' );
    delete [] list;

    for ( int ifld=0; ifld < (int)fldlist.size(); ifld++ ) {
      int32 ndims;
      double min = 0.0;
      double max = 0.0;
      DataType btype;
      DataField df;
      df.swathname = swlist[i];
      df.fieldname = fldlist[ifld];
      df.type      = type[ifld];

      // if( df.type == NOT_DEFINED )
      //    continue;

      dims = new int32[ rank[ifld] ];
      list = new char[300];
      memset( list, 0, 300 );
      nret = SWfieldinfo( swid, (char*)fldlist[ifld].c_str(),
                          &ndims, dims, &df.type, list );

      if ( ndims != 2 || dims[0] < 50 || dims[1] < 50 ) {
        delete [] dims;
        delete [] list;
        continue;
      }

      df.dimnames = splitList( list, ',' );
      delete [] list;

      for ( int idim=0; idim < ndims; idim++ )
        df.dims.push_back( dims[ndims-idim-1] );
      delete [] dims;

      switch ( df.type ) {
        case DFNT_INT8:
          min = -127;
          max = 128;
          btype = INT8;
          break;

        case DFNT_UINT8:
          min = 0;
          max = 255;
          btype = UINT8;
          break;

        case DFNT_INT16:
          {
            short lwr=32767, upr=-32767;

            // original
            //min = lwr;
            //max = upr;
            //findrange( swid, df.fieldname, lwr, upr, df.dims );

            // dls did this to fix problem reading INT16 HDF data
            findrange( swid, df.fieldname, lwr, upr, df.dims );
            min = lwr;
            max = upr;

            btype = INT16;

          }
          break;

        case DFNT_UINT16:
          {
            unsigned short lwr=65535, upr = 0;

            findrange( swid, df.fieldname, lwr, upr, df.dims );
            min = lwr;
            max = upr;

            btype = UINT16;
          }
          break;

        default:
          btype = NOT_DEFINED;
      }


      BandInf b;
      b.name   = (char*)(df.swathname+"-"+df.fieldname).c_str();
      b.ndims  = ndims;
      b.width  = df.dims[0];
      b.height = df.dims[1];
      b.max    = max;
      b.min    = min;
      b.type   = btype;

      if ( ! (b.type == NOT_DEFINED) ) {
        mBandSet.push_back( b );
        mDFList.push_back( df );
      }


    } // END DATA FIELD

    delete [] rank;
    delete [] type;

    SWdetach( swid );
  } // END SWATH

  if ( mDFList.size() == 0 )
    return false;

  // FILL IN IMAGE INFO
  if (mImageInf.mAquiDateTime.size() == 0) {
    mImageInf.mAquiDateTime      = getAquiDateTime( fname );
  }
  if (mImageInf.mCloudPct == 0) {
    mImageInf.mCloudPct          = atof( getCloudPct( fname ).c_str() );
  }
  if (mImageInf.mImageAzim == 0) {
    mImageInf.mImageAzim         = atof( getImageAzim( fname ).c_str() );
  }
  if (mImageInf.mOrigID.size() == 0) {
    mImageInf.mOrigID            = getOrigID( fname );
  }
  if (mImageInf.mInstrument_zenith == 0) {
    mImageInf.mInstrument_zenith = atof( getInstZenith( fname ).c_str() );
  }
  if (mImageInf.mSunAzim == 0) {
    mImageInf.mSunAzim           = atof( getSunAzim( fname ).c_str() );
  }
  if (mImageInf.mSunElev == 0) {
    mImageInf.mSunElev           = atof( getSunElev( fname ).c_str() );
  }
  if (mImageInf.mProjection.size() == 0) {
    mImageInf.mProjection        = getProjection( fname );
  }

  mImageInf.mLocURL = fname ;

  // 2007-06-14 dls - Removed image ID references and functions
  /*
  if (mImageInf.mImageID.size() <= 0) {
    mImageInf.mImageID            = getImageID( fname );
  }
  */
  if (mImageInf.mInstID.size() <= 0) {
    mImageInf.mInstID             = getInstID( fname );
  }
  if (mImageInf.mInstrument_azimuth <= 0) {
    mImageInf.mInstrument_azimuth = atof(getInstAz( fname ).c_str() ) ;
  }



  // LOAD THE GEOFIELDS ASSOCIATED WITH EACH BAND INTO THE DATAFIELD
  for ( int iband=0; iband < (int)mBandSet.size(); iband++ ) {
    getGeoFields( iband );
  }

  getBandMetaData();
  getULCorners();

  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;
  }

  return true;
}

void HdfEosFrmt::getGeoFields( int iband ) {
  int32           swid;
  int32           iLat = -1;
  int32           iLon = -1;
  int32           latoff[2];
  int32           lonoff[2];
  int32           latinc[2];
  int32           loninc[2];
  int32           ret;
  DataField       &df = mDFList[iband];

  /*
  if( df.dimnames.size() <= 1 )
      return;
  */

  swid = SWattach( mSWfid, (char*)df.swathname.c_str() );
  if ( swid == -1 ) return;

  for ( int ind=0; ind < (int)mGFList.size(); ind++ ) {
    if ( mGFList[ind].swathname != df.swathname )
      continue;

    if ( mGFList[ind].gfldname == "Latitude" )
      iLat = ind;

    if ( mGFList[ind].gfldname == "Longitude" )
      iLon = ind;
  }

  if ( iLat == -1 || iLon == -1 ) {
    SWdetach( swid );
    return;
  }

  // GET HORIZONTAL MAPPINGS
  ret = SWmapinfo( swid,
                   (char*)mGFList[iLon].dimnames[0].c_str(),
                   (char*)df.dimnames[0].c_str(),
                   &lonoff[0], &loninc[0] );

  // GET VERTICAL MAPPINGS
  ret = SWmapinfo( swid,
                   (char*)mGFList[iLon].dimnames[1].c_str(),
                   (char*)df.dimnames[1].c_str(),
                   &lonoff[1], &loninc[1] );

  // GET HORIZONTAL MAPPINGS
  ret = SWmapinfo( swid,
                   (char*)mGFList[iLat].dimnames[0].c_str(),
                   (char*)df.dimnames[0].c_str(),
                   &latoff[0], &latinc[0] );

  // GET VERTICAL MAPPINGS
  ret = SWmapinfo( swid,
                   (char*)mGFList[iLat].dimnames[1].c_str(),
                   (char*)df.dimnames[1].c_str(),
                   &latoff[1], &latinc[1] );

  df.geofldlat = &mGFList[iLat];
  df.geofldlon = &mGFList[iLon];

  if ( latoff[0] < 0 || latoff[0] > df.dims[0] )
    latoff[0] = 0;
  if ( latoff[1] < 0 || latoff[1] > df.dims[1] )
    latoff[1] = 0;
  if ( lonoff[0] < 0 || lonoff[0] > df.dims[0] )
    lonoff[0] = 0;
  if ( lonoff[1] < 0 || lonoff[1] > df.dims[1] )
    lonoff[1] = 0;

  df.latoff.push_back( latoff[0] );
  df.latoff.push_back( latoff[1] );
  df.lonoff.push_back( lonoff[0] );
  df.lonoff.push_back( lonoff[1] );

  df.latinc.push_back( latinc[0] );
  df.latinc.push_back( latinc[1] );
  df.loninc.push_back( loninc[0] );
  df.loninc.push_back( loninc[1] );

  SWdetach( swid );
}

void HdfEosFrmt::getLL( int iband, double &i, double &j, bool toWGS84 ) {
  double      lat;
  double      lon;
  DataField   &df = mDFList[iband];

  if ( iband < 0 || iband > (int)mBandSet.size() ||
       i < 0 || i > mBandSet[iband].width ||
       j < 0 || j > mBandSet[iband].height ) {
    i = 0;
    j = 0;
    return;
  }

  switch ( df.geofldlon->type ) {
    case DFNT_FLOAT32:
      bilinterp( (float*)df.geofldlon->geomap,
                 df.geofldlon->dims,
                 df.lonoff,
                 df.loninc,
                 (int)i,
                 (int)j,
                 lon );
      break;

    case DFNT_FLOAT64:
      bilinterp( (double*)df.geofldlon->geomap,
                 df.geofldlon->dims,
                 df.lonoff,
                 df.loninc,
                 (int)i,
                 (int)j,
                 lon );
      break;

    default:
      // CLEAN UP HERE
      break;
  }

  switch ( df.geofldlat->type ) {
    case DFNT_FLOAT32:
      bilinterp( (float*)df.geofldlat->geomap,
                 df.geofldlat->dims,
                 df.latoff,
                 df.latinc,
                 (int)i,
                 (int)j,
                 lat );
      break;

    case DFNT_FLOAT64:
      bilinterp( (double*)df.geofldlat->geomap,
                 df.geofldlat->dims,
                 df.latoff,
                 df.latinc,
                 (int)i,
                 (int)j,
                 lat );
      break;
  }


  /* CONVERT GEODEDIC TO GEOCENTRIC LAT */
  // double wgs84_C = 1.0067395;
  if ( toWGS84 ) {
    double wgs84_C = 1.00673797751;
    lat = (180.0/M_PI) * atan( ( tan(lat*(M_PI/180.0) ) * wgs84_C ) );
  }
  if (lon < 0.0) {
    lon += 360.0 ;
  }
  i = lon;
  j = lat;
}

void HdfEosFrmt::getXY( int iband, double &i, double &j ) {
  DataField   &df     = mDFList[iband];
  GeoField    &gflat  = *df.geofldlat;
  GeoField    &gflon  = *df.geofldlon;
  int         x, y;
  int         mapw, maph;
  bool        latfound, lonfound;
  double      lon, lat;
  int         imgw, imgh;

  double wgs84_C = 1.00673797751;
  std::vector<int*>     cross_set;

  if ( iband < 0 || iband > (int)mBandSet.size() ) {
    return;
  }

  if ( mBandSet[iband].proj != "" ) {
    llToXY( iband, i, j );
    return;
  }

  imgw = mBandSet[iband].width;
  imgh = mBandSet[iband].height;
  lon = i;
  lat = j;
  i = -999999.0;
  j = -999999.0;

  // CONVERT TO GEOCENTRIC COORDINATES AS THE LAT/LON GRID IS IN
  lat = atan( tan(lat / (180.0/M_PI)) / wgs84_C ) / (M_PI/180.0);

  // GET THE GEOGRAPHICAL COORDINATE MAPS FOR LAT/LON OF THIS BAND
  mapw = gflon.dims[1];
  maph = gflon.dims[0];

  // SEARCHING THE LONGITUDE TABLE HERE FOR INTERSECTION OF THE INPUT
  // LONGITUDE.  ONLY TRYING TO FIND ONE TO VERIFY IT'S ON THE GRID.
  lonfound = false;
  latfound = false;
  for ( y=0; y < maph && !lonfound; y++ ) {
    for ( x=0; x < mapw && !lonfound; x++ ) {
      // SEARCH HORIZONTAL LINE
      if ( x < mapw-1 ) {
        if ( gflon.type == DFNT_FLOAT32 ) {
          float *data = (float*)gflon.geomap;
          if ( between( data[y*mapw + x], data[y*mapw + x+1], lon ) ) {
            lonfound = true;
            continue;
          }
        } else if ( gflon.type == DFNT_FLOAT64 ) {
          double *data = (double*)gflon.geomap;
          if ( between( data[y*mapw + x], data[y*mapw + x+1], lon ) ) {
            lonfound = true;
            continue;
          }
        }
      }

      // SEARCH VERTICAL LINE
      if ( y < maph-1 ) {
        if ( gflon.type == DFNT_FLOAT32 ) {
          float *data = (float*)gflon.geomap;
          if ( between( data[y*mapw + x], data[(y+1)*mapw + x], lon ) ) {
            lonfound = true;
            continue;
          }
        } else if ( gflon.type == DFNT_FLOAT64 ) {
          double *data = (double*)gflon.geomap;
          if ( between( data[y*mapw + x], data[(y+1)*mapw + x], lon ) ) {
            lonfound = true;
            continue;
          }
        }
      }
    }
  }

  mapw = gflat.dims[1];
  maph = gflat.dims[0];
  for ( y=0; y < maph && !latfound; y++ ) {
    for ( x=0; x < mapw && !latfound; x++ ) {
      // SEARCH HORIZONTAL LINE
      if ( x < mapw-1 ) {
        if ( gflon.type == DFNT_FLOAT32 ) {
          float *data = (float*)gflat.geomap;
          if ( between( data[y*mapw + x], data[y*mapw + x+1], lat ) ) {
            latfound = true;
            continue;
          }
        } else if ( gflon.type == DFNT_FLOAT64 ) {
          double *data = (double*)gflat.geomap;
          if ( between( data[y*mapw + x], data[y*mapw + x+1], lat ) ) {
            latfound = true;
            continue;
          }
        }
      }

      // SEARCH VERTICAL LINE
      if ( y < maph-1 ) {
        if ( gflon.type == DFNT_FLOAT32 ) {
          float *data = (float*)gflat.geomap;
          if ( between( data[y*mapw + x], data[(y+1)*mapw + x], lat ) ) {
            latfound = true;
            continue;
          }
        } else if ( gflon.type == DFNT_FLOAT64 ) {
          double *data = (double*)gflat.geomap;
          if ( between( data[y*mapw + x], data[(y+1)*mapw + x], lat ) ) {
            latfound = true;
            continue;
          }
        }
      }

    }
  }

  if ( !lonfound || !latfound ) {
    return;
  }

  // SEARCH THE LAT GRID STORING ALL INTERSECTIONS.  IF ONE NOT FOUND THEN
  // INPUT COORDS NOT ON IMAGE.
  //
  // COULD THIS JUST DO 24 ITERATIONS, 12 IN X AND 12 IN Y, CHECKING FOR
  // AN INTERSECTION ON THE LINE FOR THE LENGTH OF THE GRID???
  //
  bool found = false;
  mapw = gflat.dims[1];
  maph = gflat.dims[0];

  //int track, p1, p2;
  //int *set;

  for ( y=0; y < maph; y++ ) {
    for ( x=0; x < mapw; x++ ) {
      int track, p1, p2;
      int *set;
      /*
      // SEARCH HORIZONTAL LINE
      if( x < mapw-1 )
      {
          if( gflat.type == DFNT_FLOAT32 )
          {
              float *data = (float*)gflat.geomap;
              if( between( data[y*mapw + x], data[y*mapw + x+1], lat ) )
              {
                  found = true;
                  track = df.latoff[1] + (y * df.latinc[1]);
                  p1    = df.latoff[0] + (x * df.latinc[0]);
                  p2    = df.latoff[0] + ((x+1) * df.latinc[0]);
                  getIntersect( track, p1, p2, iband, lat, false, true );
                  set = new int[2];
                  set[0] = p1;
                  set[1] = track;
                  cross_set.push_back( set );
              }
          } else if( gflat.type == DFNT_FLOAT64 ) {
              double *data = (double*)gflat.geomap;
              if( between( data[y*mapw + x], data[y*mapw + x+1], lat ) )
              {
                  found = true;
                  track = df.latoff[1] + (y * df.latinc[1]);
                  p1    = df.latoff[0] + (x * df.latinc[0]);
                  p2    = df.latoff[0] + ((x+1) * df.latinc[0]);
                  getIntersect( track, p1, p2, iband, lat, false, true );
                  set = new int[2];
                  set[0] = p1;
                  set[1] = track;
                  cross_set.push_back( set );
              }

          }
      }
      */
      // SEARCH VERTICAL LINE
      if ( y < maph-1 ) {
        if ( gflat.type == DFNT_FLOAT32 ) {
          float *data = (float*)gflat.geomap;
          if ( between( data[y*mapw + x], data[(y+1)*mapw + x], lat ) ) {
            found = true;
            track = df.latoff[0] + (x * df.latinc[0]);
            p1    = df.latoff[1] + (y * df.latinc[1]);
            p2    = df.latoff[1] + ((y+1) * df.latinc[1]);
            getIntersect( track, p1, p2, iband, lat, false, false );
            set = new int[2];
            set[0] = track;
            set[1] = p1;
            cross_set.push_back( set );
          }
        } else if ( gflat.type == DFNT_FLOAT64 ) {
          double *data = (double*)gflat.geomap;
          if ( between( data[y*mapw + x], data[(y+1)*mapw + x], lat ) ) {
            found = true;
            track = df.latoff[0] + (x * df.latinc[0]);
            p1    = df.latoff[1] + (y * df.latinc[1]);
            p2    = df.latoff[1] + ((y+1) * df.latinc[1]);
            getIntersect( track, p1, p2, iband, lat, false, false );
            set = new int[2];
            set[0] = track;
            set[1] = p1;
            cross_set.push_back( set );
          }
        }
      }
    }
  }

  if ( (int)cross_set.size() == 0 ) {
    return;
  }

  // FIND BEST INTERSECT POINT
  double testlon, testlat;
  double prevlon;
  x = cross_set[0][0];
  y = cross_set[0][1];
  prevlon = x;
  testlat = y;
  getLL( iband, prevlon, testlat );

  for ( int ipnt=1; ipnt < (int)cross_set.size(); ipnt++ ) {
    testlon = cross_set[ipnt][0];
    testlat = cross_set[ipnt][1];
    getLL( iband, testlon, testlat );
    if ( fabs( testlon - lon ) < fabs( prevlon - lon ) ) {
      prevlon = testlon;
      x = cross_set[ipnt][0];
      y = cross_set[ipnt][1];
    }
  }

  trackLineTo( x, y, iband, lat, lon );

  for ( int ivec=0; ivec < (int)cross_set.size(); ivec++ )
    delete cross_set[ivec];

  i = x;
  j = y;
}

void HdfEosFrmt::getIntersect( int &track,
                               int &p1,
                               int p2,
                               int iband,
                               double searchVal,
                               bool,
                               bool horizontal ) {
  int l, m, r;
  double lon, lat;
  double val;
  bool ascending = true;

  if ( p2 < p1 ) {
    int tmp = p1;
    p1 = p2;
    p2 = tmp;
  }

  l = p1;
  r = p2;
  m = (r - l) / 2;

  // CHECK FOR TRACK VALUE DIRECTION
  if ( horizontal ) {
    double lon1, lat1, lon2, lat2;
    lon1 = l;
    lat1 = lat2 = track;
    lon2 = r;
    getLL( iband, lon1, lat1, false );
    getLL( iband, lon2, lat2, false );
    ascending = ( lat1 < lat2 );
  } else {
    double lon1, lat1, lon2, lat2;
    lat1 = l;
    lon1 = lon2 = track;
    lat2 = r;
    getLL( iband, lon1, lat1, false );
    getLL( iband, lon2, lat2, false );
    ascending = ( lat1 < lat2 );
  }

  while ( m > 0 ) {
    if ( horizontal ) {
      lon = l + m;
      lat = track;
    } else {
      lon = track;
      lat = l + m;
    }

    getLL( iband, lon, lat, false );

    val = lat;
    if ( ( ascending && val < searchVal ) ||
         ( !ascending && val > searchVal ) )
      l += m;
    else
      r -= m;

    m = (r - l) / 2;
  }

  p1 = l;
  if ( l != r ) {
    double lon1, lat1, lon2, lat2;
    double val2;
    if ( horizontal ) {
      lon1 = l;
      lat1 = track;
      lon2 = r;
      lat2 = track;
    } else {
      lon1 = track;
      lat1 = l;
      lon2 = track;
      lat2 = r;
    }

    getLL( iband, lon1, lat1, false );
    getLL( iband, lon2, lat2, false );

    val = lat1;
    val2 = lat2;
    if ( fabs( val - searchVal ) > fabs( val2 - searchVal ) )
      p1 = r;
  }

  return;
}

void HdfEosFrmt::trackLineTo( int &x,
                              int &y,
                              int iband,
                              double trackVal,
                              double searchVal ) {
  int a = 0, b = 0;
  int x1 = 0, y1 = 0, x2 = 0, y2 = 0;
  bool found = false;
  int cnt = 0;

  double bestlat1 = 0, bestlat2 = 0;
  double   testlon[3][3] = {
    { -999.0, -999.0, -999.0},
    { -999.0, -999.0, -999.0},
    { -999.0, -999.0, -999.0}
  };

  double   testlat[3][3] = {
    { -999.0, -999.0, -999.0},
    { -999.0, -999.0, -999.0},
    { -999.0, -999.0, -999.0}
  };

  while ( !found ) {
    for ( a = 0; a < 3; a++ ) {
      for ( b = 0; b < 3; b++ ) {
        testlat[a][b] = y + a-1;
        testlon[a][b] = x + b-1;
        getLL( iband, testlon[a][b], testlat[a][b], false );
      }
    }

    double londiff = fabs( testlon[1][1] - searchVal );

    found = true;
    for ( a = 0; found && a < 3; a++ ) {
      for ( b = 0; b < 3; b++ ) {
        if ( a == 1 && b == 1 )
          continue;

        if ( fabs( testlon[a][b] - searchVal ) < londiff ) {
          found = false;
          continue;
        }
      }
    }

    bestlat1 = bestlat2 = 999.0;
    for ( a=0; a < 3; a++ ) {
      for ( b=0; b < 3; b++ ) {
        if ( a == 1 && b == 1 )
          continue;

        double latdiff = fabs( testlat[a][b] - trackVal );
        if ( latdiff < bestlat1 ) {
          bestlat2 = bestlat1;
          x2 = x1;
          y2 = y1;
          x1 = b;
          y1 = a;
          bestlat1 = latdiff;
        } else if ( latdiff < bestlat2 ) {
          x2 = b;
          y2 = a;
          bestlat2 = latdiff;
        }
      }
    }



    if ( fabs( searchVal - testlon[y1][x1] ) >
         fabs( searchVal - testlon[y2][x2] ) ) {
      x1 = x2;
      y1 = y2;
    }

    /*
    if( fabs( searchVal - testlon[1][1] ) <
        fabs( searchVal - testlon[y1][x1] ) )
    {
        found = true;
    } else {
        x += x1-1;
        y += y1-1;
    }
    */

    x += x1-1;
    y += y1-1;

    cnt++;
  }
}

void * HdfEosFrmt::getBand( int iband,
                            int x,
                            int y,
                            int w,
                            int h,
                            int w2,
                            int h2,
                            ScaleType st ) const {
  unsigned char* data = NULL;

  if ( w <= 0 )
    w = mBandSet[iband].width;
  if ( h <= 0 )
    h = mBandSet[iband].height;

  if ( w2 == 0 ) return data; // DIVIDE BY ZERO

  int32 skip = w / w2;
  if ( skip == 0 )
    skip = 1;

  w = (w / skip) * skip;
  h = (h / skip) * skip;

  int nw = w / skip;
  int nh = h / skip;

  int32 start[]     = { y, x};
  int32 scales[]    = { skip, skip};
  int32 edge[]      = { nh, nw};
  int32 swid;
  int     dSize;

  dSize = ImageFormat::dtypeSize( mBandSet[iband].type );
  data = new unsigned char[(nw * nh) * dSize];

  std::string swname  = mDFList[iband].swathname;
  std::string fldname = mDFList[iband].fieldname;

  swid = SWattach( mSWfid, (char*)swname.c_str() );
  if ( swid == -1 ) {
    delete [] data;
    return NULL;
  }

  SWreadfield( swid, (char*)fldname.c_str(), start, scales, edge, data );
  SWdetach( swid );

  if ( nw != w2 || nh != h2 )
    scaleData( (char**)&data, nw, nh, w2, h2, mBandSet[iband].type, st );
  return data;
}

void HdfEosFrmt::parseGeoFlds( int32 swid, std::string swname ) {
  int32 nret, ngeofld, bufsz;
  int32 *rank, *type;
  char  *list;
  std::vector<std::string> gflist;

  ngeofld = SWnentries( swid, HDFE_NENTGFLD, &bufsz );

  list = new char[bufsz];
  rank = new int32[bufsz];
  type = new int32[bufsz];
  nret = SWinqgeofields( swid, list, rank, type );
  if ( nret < 1 )
    return;

  gflist = splitList( list, ',' );
  delete [] list;

  for ( int igfld=0; igfld < (int)gflist.size(); igfld++ ) {
    int32           *dims;
    int32           ndims;
    int             datasz = 0;
    GeoField        gf;

    gf.swathname = swname;
    gf.gfldname  = gflist[igfld];

    dims = new int32[ rank[igfld] ];
    list = new char[1024];
    memset( list, 0, 1024 );
    nret = SWfieldinfo( swid, (char*)gflist[igfld].c_str(),
                        &ndims, dims, &gf.type, list );

    if ( nret == -1 ) {
      delete [] list;
      delete [] dims;
      continue;
    }

    gf.dimnames = splitList( list, ',' );
    delete [] list;

    for ( int idim=0; idim < ndims; idim++ )
      gf.dims.push_back( dims[idim] );
    delete [] dims;


    switch ( gf.type ) {
      case DFNT_FLOAT32:
        datasz = 4;
        break;

      case DFNT_FLOAT64:
        datasz = 8;
        break;
    }

    if ( datasz == 0 )
      continue;

    gf.geomap = new unsigned char[ gf.dims[0] * gf.dims[1] * datasz ];

    int32 start[] = { 0, 0};
    int32 scale[] = { 1, 1};
    int32 edge[]  = { gf.dims[0], gf.dims[1]};

    SWreadfield( swid, (char*)gflist[igfld].c_str(), start, scale, edge, gf.geomap );
    mGFList.push_back( gf );
  }

  delete [] rank;
  delete [] type;
}

template<class DType>
void HdfEosFrmt::findrange( int32 swid,
                            std::string field,
                            DType &lwr,
                            DType &upr,
                            std::vector<int32> dims ) {
  int32 ret;

  if ( (int)dims.size() < 1 )
    return;

  int32 start[2]  = { 0, 0};
  int32 scales[2] = { 1, 1};
  int32 edge[2]   = { dims[1], dims[0]};

  int32 bufsz = dims[1] * dims[0];

  DType* buf = new DType[ bufsz ];
  if ( buf == NULL ) {
    char err[256];
    sprintf( err, "buf is NULL %d", (int)bufsz );
    QMessageBox::information( NULL, "", err );
    return;
  }
  ret = SWreadfield( swid, (char*)field.c_str(), start, scales, edge, (unsigned char*)buf );
  maxmin( buf, bufsz, lwr, upr );
  delete [] buf;

}


std::string HdfEosFrmt::getAquiDateTime( const std::string &filename ) const {
  std::string strtod;
  std::string strcal;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( filename, "coremetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "TIMEOFDAY" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      strtod = elem->getValue();
  }

  tree = p.getObject( "CALENDARDATE" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      strcal = elem->getValue();
  }

  // check for invalid characters in TIMEOFDAY field
  bool valid = true ;
  char tmpstr[2] = " " ;
  unsigned int i ;
  for (i = 0 ; i < strtod.length() ; i++) {
    tmpstr[0] = strtod[i] ;
    if (!(tmpstr[0] == '\"' ||
          tmpstr[0] == '.'  ||
          tmpstr[0] == 'Z'  ||
          (tmpstr[0] >= '0' && tmpstr[0] <= '9') ) ) {
      valid = false ;
    }
  }
  // If there are no invalid chars in the TIMEOFDAY field
  // reformat the string
  if (valid) {
    std::string h,m,s;
    if ( strtod.length() > 0 && strtod[0] == '\"' )
      strtod = strtod.substr( 1, strtod.length()-2 );
    h = strtod.substr( 0, 2 );
    m = strtod.substr( 2, 2 );
    s = strtod.substr( 4, 2 );

    strtod = h + ":" + m + ":" + s + "-00:00";
  }


  // check for invalid characters in CALENDARDATE field
  for (i = 0 ; i < strcal.length() ; i++) {
    tmpstr[0] = strcal[i] ;
    if (!(tmpstr[0] == '\"' ||
          tmpstr[0] == '.'  ||
          (tmpstr[0] >= '0' && tmpstr[0] <= '9') ) ) {
      valid = false ;
    }
  }
  // If there are no invalid chars in the CALENDARDATE field
  // reformat the string
  if (valid) {
    if ( strcal.length() > 0 && strcal[0] == '\"' ) {
      strcal = strcal.substr( 1, strcal.length()-2 );
    }

    strcal = strcal.substr( 0, 4 ) + "-" +
             strcal.substr( 4, 2 ) + "-" +
             strcal.substr( 6, 2 ) ;

  }

  if (valid) {
    // return ISO 1806 compliant
    return( strcal + "T" + strtod );
  }
  else {
    // return exactly what was stored in the metadata
    return( strcal + " " + strtod );
  }
}

std::string HdfEosFrmt::getOrigID( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "IDOFASTERGDSDATAGRANULE" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }

  if ( val.length() > 0 && val[0] == '\"' )
    val = val.substr( 1, val.length()-2 );

  return val;
}

std::string HdfEosFrmt::getInstZenith( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "POINTINGANGLE" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 2 );
    if ( elem )
      val = elem->getValue();
  }

  return val;
}

std::string HdfEosFrmt::getCloudPct( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "SCENECLOUDCOVERAGE" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }

  return val;
}

std::string HdfEosFrmt::getSunAzim( const std::string &fname ) const {
  std::string val;
  int i1;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "SOLARDIRECTION" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }

  i1 = val.find( '(' )+1;
  val = val.substr( i1, val.find( ',' )-i1 );
  return val;
}

std::string HdfEosFrmt::getSunElev( const std::string &fname ) const {
  std::string val;
  int i1;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "SOLARDIRECTION" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }

  i1 = val.find( ',' )+1;
  val = val.substr( i1, val.find( ')' )-i1 );
  return val;
}

std::string HdfEosFrmt::getImageAzim( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "MAPORIENTATIONANGLE" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }

  return val;
}

std::string HdfEosFrmt::getProjection( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "coremetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "MAPPROJECTIONNAME" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }

  if ( val.length() > 0 && val[0] == '\"' )
    val = val.substr( 1, val.length()-2 );

  return val;
}



// 2007-06-14 dls - Removed image ID references and functions
/*
std::string   HdfEosFrmt::getImageID( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "productmetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "ASTERSCENEID" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem ) {
      val = elem->getValue();
      if ( val.length() > 0 && val[0] == '(' ) {
        val = val.substr( 1, val.length()-2 );
        unsigned int i = 0 ;
        QString build = "Path=" ;
        while (i < val.length() && val[i++] != ',') {
          if (val[i] >= '0' && val[i] <= '9') {
            build += val[i] ;
          }
        }
        build += " Row=" ;
        while (i < val.length() && val[++i] != ',') {
          if (val[i] >= '0' && val[i] <= '9') {
            build += val[i] ;
          }
        }
        build += " View=" ;
        while (i < val.length() && val[++i] != ',') {
          if (val[i] >= '0' && val[i] <= '9') {
            build += val[i] ;
          }
        }
        val = build.ascii() ;

      }
    }
  }

  return val;
}
*/

std::string   HdfEosFrmt::getInstID( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "coremetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "INSTRUMENTSHORTNAME" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 0 );
    if ( elem )
      val = elem->getValue();
  }

  if ( val.length() > 0 && val[0] == '\"' )
    val = val.substr( 1, val.length()-2 );

  return val;
}

std::string   HdfEosFrmt::getInstAz( const std::string &fname ) const {
  std::string val;
  const ODLTree *tree;

  std::string data = getHdfGblAttr( fname, "coremetadata.0" );
  if ( data == "" )
    return "";

  ODLParser p( data );
  tree = p.getObject( "SCENEORIENTATIONANGLE" );
  if ( tree ) {
    ODLElement *elem = tree->getElem( 1 );
    if ( elem )
      val = elem->getValue();
  }
  return val;
}





std::string HdfEosFrmt::getHdfGblAttr(
                                     const std::string &fname,
                                     const std::string &attrname ) const {
  int32           sd_id;
  //int32           sds_id;
  int32           status;
  int32           dtype;                      // DATA TYPE
  int32           count;                      // NUMBER OF ENTRIES
  int32           iattr;                      // ATTRIBUTE INDEX
  int32           ndset;                      // # OF DATASETS
  int32           nattr;                      // # OF ATTRIBUTES
  char            attrnamebuf[MAX_NC_NAME];   // ATTRIBUTE NAME
  int8            *attrbuf;                   // ATTRIBUTE TEXT BUFFER

  // Open the HDF file
  sd_id   = SDstart( fname.c_str(), DFACC_READ );
  if ( !sd_id ) return "";

  status  = SDfileinfo( sd_id, &ndset, &nattr );
  iattr   = SDfindattr( sd_id, attrname.c_str() );
  if ( iattr < 0 ) return "";

  status  = SDattrinfo( sd_id, iattr, attrnamebuf, &dtype, &count );
  attrbuf = new int8[ count * DFKNTsize( dtype ) ];

  status  = SDreadattr( sd_id, iattr, attrbuf );

  // Dispose of the file identifier to close the file.
  status = SDend( sd_id );

  std::string ret = attrbuf;
  delete [] attrbuf;
  return ret;
}

// PRETTY ASTER SPECIFIC HERE
void HdfEosFrmt::getBandMetaData( ) const {
  std::string         IMGD = "ImageData";
  std::string         bid;
  std::string         swath;
  std::string         swathFull;
  int                 pos;
  std::string         val;
  const ODLTree       *tree;
  std::string         data;
  double              rot = 0.0;
  double              ulx = 0.0, uly = 0.0;

  std::string         prod_v;
  std::string         prod_s;
  std::string         prod_t;
  std::string         prod_0;

  ODLParser           *pv;
  ODLParser           *ps;
  ODLParser           *pt;
  ODLParser           *p0;

  prod_v = getHdfGblAttr( mFName, std::string( "productmetadata.v" ) );
  prod_s = getHdfGblAttr( mFName, std::string( "productmetadata.s" ) );
  prod_t = getHdfGblAttr( mFName, std::string( "productmetadata.t" ) );
  prod_0 = getHdfGblAttr( mFName, std::string( "productmetadata.0" ) );

  pv = new ODLParser( prod_v );
  ps = new ODLParser( prod_s );
  pt = new ODLParser( prod_t );
  p0 = new ODLParser( prod_0 );

  if ( p0 ) {
    tree = p0->getObject( std::string( "MAPORIENTATIONANGLE" ) );
    if ( tree ) {
      ODLElement *elem = tree->getElem( 1 );
      if ( elem )
        rot = atof( elem->getValue().c_str() );
    } else {
      tree = p0->getObject( std::string( "SCENEORIENTATIONANGLE" ) );
      if ( tree ) {
        ODLElement *elem = tree->getElem( 1 );
        if ( elem )
          rot = 0-atof( elem->getValue().c_str() );
      }
    }

    tree = p0->getObject( std::string( "SCENEFOURCORNERS" ) );
    if ( tree ) {
      const std::vector<ODLTree> &children = tree->getChildSet();
      for ( int ichld=0; ichld < (int)children.size(); ichld++ ) {
        if ( children[ichld].getName() == "UPPERLEFT" ) {
          ODLElement *elem = children[ichld].getElem( 1 );
          if ( elem ) {
            int l, r;
            std::string val = elem->getValue();
            l = val.find( "(" ) + 1;
            r = val.find( ")" );
            val = val.substr( l, r-l );
            l = 0;
            r = val.find( "," );
            uly = atof( val.substr( l, r-l ).c_str() );
            val = val.substr( r+1, val.length()-r );
            ulx = atof( val.c_str() );
          }
        }
      }
      /*
      ODLElement *elem = tree->getElem( 1 );
      if( elem )
          rot = atof( elem->getValue().c_str() );
      */
    }
  }

  for ( int iband=0; iband < (int)mBandSet.size(); iband++ ) {
    BandInf &binf = (BandInf&)mBandSet[iband];
    const ODLTree *mptree, *uztree;

    binf.ulx = ulx;
    binf.uly = uly;

    pos = binf.name.find_last_of( IMGD ) + 1;
    if ( pos <= 0 ) continue;
    bid = binf.name.substr( pos );
    if ( bid.length() < 1 ) continue;

    if ( (int)binf.name.find( "VNIR" ) != -1 ) {
      mptree = pv->getObject( std::string( "MPMETHOD" ) + bid );
      uztree = pv->getObject( std::string( "UTMZONECODE" ) + bid );
      binf.res = 15;
    } else if ( (int)binf.name.find( "SWIR" ) != -1 ) {
      mptree = ps->getObject( std::string( "MPMETHOD" ) + bid );
      uztree = ps->getObject( std::string( "UTMZONECODE" ) + bid );
      binf.res = 30;
    } else if ( (int)binf.name.find( "TIR" ) != -1 ) {
      mptree = pt->getObject( std::string( "MPMETHOD" ) + bid );
      uztree = pt->getObject( std::string( "UTMZONECODE" ) + bid );
      binf.res = 90;
    } else {
      continue;
    }

    if ( mptree ) {
      ODLElement *elem = mptree->getElem( 1 );
      if ( elem ) {
        std::string val = elem->getValue();
        if ( (int)val.find( "UTM" ) >= 0 )
          binf.proj = "utm";
        else if ( (int)val.find( "EQRECT" ) >= 0 )
          binf.proj = "";
        else if ( (int)val.find( "LAMCC" ) >= 0 )
          binf.proj = "";
        else if ( (int)val.find( "SOM" ) >= 0 )
          binf.proj = "";
        else if ( (int)val.find( "PS" ) >= 0 )
          binf.proj = "";
        else
          binf.proj = "";
      }
    }

    if ( uztree && binf.proj != "" ) {
      ODLElement *elem = uztree->getElem( 1 );
      if ( elem )
        binf.UTMZ = atoi( elem->getValue().c_str() );
    }

    binf.rot = rot;
    binf.parms.push_back( "ellps=WGS84" );
  }

  delete pv;
  delete ps;
  delete pt;
  delete p0;
}

template<class DType>
void HdfEosFrmt::bilinterp( DType  *map,
                            std::vector<int32>  dims,
                            std::vector<int32>  offset,
                            std::vector<int32>  inc,
                            int    x,
                            int    y,
                            double &res ) {
  if ( dims[0] < 1 || dims[1] < 1 )
    return;

  int minX = offset[0];
  int minY = offset[1];

  int32 xinc = inc[1];
  int32 yinc = inc[0];
  int32 xdim = dims[1];
  int32 ydim = dims[0];

  int maxX = xdim*xinc + minX;
  int maxY = ydim*yinc + minY;

  if ( x < minX || x > maxX ||
       y < minY || y > maxY )
    return;

  int gridX = (int)(((float)(x-minX) / (float)(maxX-minX)) * xdim);
  int gridY = (int)(((float)(y-minY) / (float)(maxY-minY)) * ydim);

  int mapY = gridY * xdim;
  int mapYp1 = (gridY+1) * xdim;

  DType dnA = map[mapY   + gridX];
  DType dnB = map[mapY   + gridX+1];
  DType dnC = map[mapYp1 + gridX];
  DType dnD = map[mapYp1 + gridX+1];

  double deltaX = (x - (gridX * xinc + offset[0])) / (double)xinc;
  double deltaX2 = 1.0-deltaX;
  double deltaY = (y - (gridY * yinc + offset[1])) / (double)yinc;

  res = (deltaX*dnB + deltaX2*dnA) * (1-deltaY) +
        (deltaX*dnD + deltaX2*dnC) * deltaY;
}

---

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