Raster.icc

//-*- C++ -*-
//-*- coding: utf-8 -*-
//
// Author: Marco Lavalle
// Original code: Joshua Cohen
// Copyright 2018
//
 
#if !defined(ISCE_IO_RASTER_ICC)
#error "Raster.icc is an implementation detail of class Raster"
#endif
 
#include <iostream>
#include <isce3/except/Error.h>
 
//Constructor to map a raw stream to a Raster
inline void isce3::io::Raster::initFromPointer(void* ptr,  //Raw pointer
                        GDALDataType dtype, //GDAL datatype
                        size_t width,  //width of raster
                        size_t length, //length of raster
                        size_t pixeloffset, //pixeloffset in bytes
                        size_t lineoffset) //lineoffset in bytes
{
    //Get mem driver
    GDALDriver* memDriver = (GDALDriver*) GDALGetDriverByName("MEM");
 
    //Create an empty dataset of requested size
    //Update with pyre error logging
    dataset(memDriver->Create("MEM::::", width, length, 0, dtype, nullptr));
 
    //Add the band information as needed
    char** memdsOptions = nullptr;
 
    {
        char szTmp[64];
 
        //Set data pointer
        memset(szTmp, 0, sizeof(szTmp));
        CPLPrintPointer(szTmp, ptr, sizeof(szTmp));
        memdsOptions = CSLSetNameValue(memdsOptions, "DATAPOINTER", szTmp);
 
        //Set line offset
        memset(szTmp, 0, sizeof(szTmp));
        CPLPrintUIntBig(szTmp, lineoffset, sizeof(szTmp));
        memdsOptions = CSLSetNameValue(memdsOptions, "LINEOFFSET", szTmp);
 
        //Set pixel offset
        memset(szTmp, 0, sizeof(szTmp));
        CPLPrintUIntBig(szTmp, pixeloffset, sizeof(szTmp));
        memdsOptions = CSLSetNameValue(memdsOptions, "PIXELOFFSET", szTmp);
    }
 
    dataset()->AddBand(dtype, memdsOptions);
}
 
 
// Constructor for a 1 band dataset from isce3::core::Matrix<T> view type
template <typename Derived>
isce3::io::Raster::Raster(Eigen::PlainObjectBase<Derived>& view)
{
    using Scalar = typename Eigen::PlainObjectBase<Derived>::value_type;
 
    //Get the packing. Update with pyre error logging.
    if (!view.IsRowMajor)
    {
        throw std::runtime_error("Input view is not packed in row major order");
    }
 
    //Size of each element
    size_t bytesperunit = sizeof(Scalar);
 
    //Pointer and offset math
    size_t pixeloffset = (char*) &view(0, 1) - (char*) &view(0, 0);
    size_t lineoffset = (char*) &view(1, 0) - (char*) &view(0, 0);
 
    //Update with pyre error logging
    if ((pixeloffset < bytesperunit) || (lineoffset < bytesperunit))
    {
        throw std::runtime_error("Something looks fishy in the offset estimation");
    }
 
    initFromPointer(view.data(),
                    asGDT<Scalar>,
                    view.cols(),
                    view.rows(),
                    pixeloffset,
                    lineoffset);
}
 
 
// Copy constructor
inline isce3::io::Raster& isce3::io::Raster::operator=(const Raster &rhs) {
 
    if (not rhs._owner) {
        throw isce3::except::RuntimeError(ISCE_SRCINFO(), "cannot copy non-owning raster");
    }
 
    if (_owner) {
        GDALClose(_dataset);
    }
 
    dataset( rhs._dataset );      // weak-copy pointer
    dataset()->Reference();       // increment GDALDataset reference counter
    return *this;
}
 
 

inline void isce3::io::Raster::open(const std::string &fname,
                                   GDALAccess access=GA_ReadOnly) {
  GDALClose( _dataset );
  dataset( static_cast<GDALDataset*>(GDALOpenShared( fname.c_str(), access )) );
}
 

inline void isce3::io::Raster::addRasterToVRT(const isce3::io::Raster& rast) {
    if (!match(rast)) // if rasters do not have the same size
        throw isce3::except::LengthError(ISCE_SRCINFO(), "Rasters have different sizes.");
 
    for (size_t b = 1; b <= rast.numBands(); ++b)       // for each band in input Raster
        addBandToVRT(rast.dataset()->GetRasterBand(b)); // add GDALRasterBand to VRT
}
 

inline void isce3::io::Raster::addBandToVRT(GDALRasterBand * inBand) {
    if (strcmp(dataset()->GetDriverName(), "VRT")) // dataset is not VRT
        throw isce3::except::RuntimeError(ISCE_SRCINFO(), "GDAL driver must be VRT.");
 
    dataset()->AddBand(inBand->GetRasterDataType(), nullptr);        // add band to VRT dataset
    GDALRasterBand *rastBand = dataset()->GetRasterBand(numBands()); // get last band
    VRTAddSimpleSource(rastBand, inBand,
                       0, 0, inBand->GetXSize(), inBand->GetYSize(),
                       0, 0, inBand->GetXSize(), inBand->GetYSize(),
                       nullptr, VRT_NODATA_UNSET);
}
 

inline void isce3::io::Raster::addRawBandToVRT(const std::string &fname, GDALDataType dtype) {
    if (strcmp(dataset()->GetDriverName(), "VRT"))
        throw isce3::except::LengthError(ISCE_SRCINFO(), "GDAL driver must be VRT.");
 
    char** papszOptions = nullptr;
    std::string srcName = fname.substr(0, fname.find_last_of("."));
    papszOptions = CSLAddNameValue(papszOptions, "subclass", "VRTRawRasterBand");
    papszOptions = CSLAddNameValue(papszOptions, "SourceFilename", srcName.c_str());
    papszOptions = CSLAddNameValue(papszOptions, "RelativeToVRT", "true");
    dataset()->AddBand(dtype, papszOptions);
    CSLDestroy(papszOptions);
}
 
 
 
 
/* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 *                                          PIXEL OPERATIONS
 * = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 */
template<typename T>
void isce3::io::Raster::getSetValue(T & buffer,           // i/o buffer of type T
                                   size_t xidx,          // 0-indexed line within band
                                   size_t yidx,          // 0-indexed column within band
                                   size_t band,          // 1-indexed band number
                                   GDALRWFlag iodir) {   // i/o direction (GF_Read or GF_Write)
 
    auto iostat = _dataset->GetRasterBand(band)->RasterIO(
            iodir, xidx, yidx, 1, 1, &buffer, 1, 1, asGDT<T>, 0, 0);
 
    if (iostat != CPLE_None) // RasterIO returned error
        std::cout << "In isce3::io::Raster::getSetValue() - error in RasterIO." << std::endl;
}
 

template<typename T>
void isce3::io::Raster::getValue(T &buffer, size_t xidx, size_t yidx, size_t band) {
    getSetValue(buffer, xidx, yidx, band, GF_Read);
}
 

template<typename T>
void isce3::io::Raster::setValue(T &buffer, size_t xidx, size_t yidx, size_t band) {
    getSetValue(buffer, xidx, yidx, band, GF_Write);
}
 
 
/* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 *                                          LINE OPERATIONS
 * = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 */
template<typename T>
void isce3::io::Raster::getSetLine(T *buffer,          // i/o buffer of type T
                                  size_t yidx,        // line location within band (0-indexed)
                                  size_t iowidth,     // width of i/o buffer
                                  size_t band,        // band number (1-indexed)
                                  GDALRWFlag iodir) { // i/o direction (GF_Read or GF_Write)
 
    size_t rdwidth = std::min(iowidth, width()); // read the requested iowidth up to width()
    auto iostat = _dataset->GetRasterBand(band)->RasterIO(iodir, 0, yidx, rdwidth, 1, buffer,
                                                          rdwidth, 1, asGDT<T>, 0, 0);
 
    if (iostat != CPLE_None) // RasterIO returned errors
        std::cout << "In isce3::io::Raster::get/setLine() - error in RasterIO." << std::endl;
}
 

template<typename T>
void isce3::io::Raster::getLine(T *buffer, size_t yidx, size_t iowidth, size_t band) {
    getSetLine(buffer, yidx, iowidth, band, GF_Read);
}
 
 
// Get a line of data for given line position and band to std::vector.
template<typename T>
void isce3::io::Raster::getLine(std::vector<T> &buffer, size_t yidx, size_t band) {
    getLine(buffer.data(), yidx, buffer.size(), band);
}
 
 
// Set a line of data for given line position and band from a std::valarray.
template<typename T>
void isce3::io::Raster::getLine(std::valarray<T> &buffer, size_t yidx, size_t band) {
    getLine(&buffer[0], yidx, buffer.size(), band);
}
 

template<typename T>
void isce3::io::Raster::setLine(T* buffer, size_t yidx, size_t iowidth, size_t band) {
    getSetLine(buffer, yidx, iowidth, band, GF_Write);
}
 
 
// Set a line of data for given line position and band from a std::vector.
template<typename T>
void isce3::io::Raster::setLine(std::vector<T> &buffer, size_t yidx, size_t band) {
    setLine(buffer.data(), yidx, buffer.size(), band);
}
 
 
// Set a line of data for given line position and band from a std::valarray.
template<typename T>
void isce3::io::Raster::setLine(std::valarray<T> &buffer, size_t yidx, size_t band) {
    setLine(&buffer[0], yidx, buffer.size(), band);
}
 
 
/* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 *                                      BLOCK OPERATIONS
 * = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 */
// Get or set a 2D data block of width iowidth and length iolength given xidx/yidx positions in band
// and to/from 1D contigous buffer of width iowidth*iolength.
template<typename T>
void isce3::io::Raster::getSetBlock(T *buffer,            // i/o buffer of size iowidth*iolength
                                   size_t xidx,          // column location within band (0-indexed)
                                   size_t yidx,          // row location within band (0-indexed)
                                   size_t iowidth,       // requested width of block of data
                                   size_t iolength,      // requested length of block of data
                                   size_t band,          // band number (1-indexed)
                                   GDALRWFlag iodir) {   // i/o direction (GF_Read or GF_Write)
 
    auto iostat = _dataset->GetRasterBand(band)->RasterIO(iodir, xidx, yidx, iowidth,
                                                          iolength, buffer, iowidth,
                                                          iolength, asGDT<T>,
                                                          0, 0);
 
    if (iostat != CPLE_None) // RasterIO returned errors
        std::cout << "In isce3::io::Raster::get/setValue() - error in RasterIO." << std::endl;
}
 

template<typename T>
void isce3::io::Raster::getBlock(T *buffer, size_t xidx, size_t yidx,
                                size_t iowidth, size_t iolength, size_t band) {
    getSetBlock(buffer, xidx, yidx, iowidth, iolength, band, GF_Read);
}
 
// Get a block of data starting from given x/y position and size in band to std::vector.
template<typename T>
void isce3::io::Raster::getBlock(std::vector<T> &buffer, size_t xidx, size_t yidx,
                                size_t iowidth, size_t iolength, size_t band) {
 
    if ((iolength * iowidth) > buffer.size()) // requested block is larger than buffer
        throw isce3::except::LengthError(ISCE_SRCINFO(), "Requested more elements than buffer size.");
 
    getBlock(buffer.data(), xidx, yidx, iowidth, iolength, band);
 
    if ((iolength * iowidth) < buffer.size())  // buffer is larger than requested block
        std::cout << "WARNING: In isce3::io::Raster::getBlock() - Requested fewer elements than buffer size. " << std::endl;
}
 
 
// Get a block of data starting from given x/y position and size in band to std::valarray.
template<typename T>
void isce3::io::Raster::getBlock(std::valarray<T> &buffer, size_t xidx, size_t yidx,
                                size_t iowidth, size_t iolength, size_t band) {
 
    if ((iolength * iowidth) > buffer.size()) // requested block is larger than buffer
        throw isce3::except::LengthError(ISCE_SRCINFO(), "Requested more elements than buffer size.");
 
    getBlock(&buffer[0], xidx, yidx, iowidth, iolength, band);
 
    if ((iolength * iowidth) < buffer.size()) // buffer is larger than requested block
        std::cout << "WARNING: In isce3::io::Raster::getBlock() - Requested fewer elements than buffer size. " << std::endl;
}
 
 
// Get a block of data starting from given x/y position in band to Matrix<T>
template<typename T>
void isce3::io::Raster::getBlock(isce3::core::Matrix<T> &mat, size_t xidx, size_t yidx,
                                size_t band) {
    getBlock(mat.data(), xidx, yidx, mat.width(), mat.length(), band);
}
 
// Get a block of data starting from given x/y position in band to EArray2D<T>
template<typename T>
void isce3::io::Raster::getBlock(isce3::core::EArray2D<T> &mat, size_t xidx, size_t yidx,
                                  size_t band) {
      getBlock(mat.data(), xidx, yidx, mat.cols(), mat.rows(), band);
}
 
// Get a block of data starting from given x/y position in band to EMatrix2D<T>
template<typename T>
void isce3::io::Raster::getBlock(isce3::core::EMatrix2D<T> &mat, size_t xidx, size_t yidx,
                                    size_t band) {
      getBlock(mat.data(), xidx, yidx, mat.cols(), mat.rows(), band);
}

template<typename T>
void isce3::io::Raster::setBlock(T *buffer, size_t xidx, size_t yidx,
                                size_t iowidth, size_t iolength, size_t band) {
    getSetBlock(buffer, xidx, yidx, iowidth, iolength, band, GF_Write);
}
 
 
// Set a block of data for given x/y position and size in band from a std::vector.
template<typename T>
void isce3::io::Raster::setBlock(std::vector<T> &buffer, size_t xidx, size_t yidx,
                                size_t iowidth, size_t iolength, size_t band) {
 
    if ((iolength * iowidth) > buffer.size()) // buffer is smaller than requested block
        throw isce3::except::LengthError(ISCE_SRCINFO(), "Requested more elements than buffer size.");
 
    setBlock(buffer.data(), xidx, yidx, iowidth, iolength, band);
 
    if ((iolength * iowidth) < buffer.size())  // buffer is larger than requested block
        std::cout << "WARNING: In isce3::io::Raster::getBlock() - Requested fewer elements than buffer size. " << std::endl;
}
 
 
// Set a block of data for given x/y position and size in band from a std::vector.
template<typename T>
void isce3::io::Raster::setBlock(std::valarray<T> &buffer, size_t xidx, size_t yidx,
                                size_t iowidth, size_t iolength, size_t band) {
 
    if ((iolength * iowidth) > buffer.size()) // buffer is smaller than requested block
        throw isce3::except::LengthError(ISCE_SRCINFO(), "Requested more elements than buffer size.");
 
    setBlock(&buffer[0], xidx, yidx, iowidth, iolength, band);
 
    if ((iolength * iowidth) < buffer.size())  // buffer is larger than requested block
        std::cout << "WARNING: In isce3::io::Raster::getBlock() - Requested fewer elements than buffer size. " << std::endl;
}
 
 
// Set a block of data for given x/y position and band from Matrix<T>
template<typename T>
void isce3::io::Raster::setBlock(isce3::core::Matrix<T>& mat, size_t xidx, size_t yidx,
                                size_t band) {
    setBlock(mat.data(), xidx, yidx, mat.width(), mat.length(), band);
}
 
// Set a block of data for given x/y position and band from EArray2D<T>
template<typename T>
void isce3::io::Raster::setBlock(isce3::core::EArray2D<T>& mat, size_t xidx, size_t yidx,
                                  size_t band) {
      setBlock(mat.data(), xidx, yidx, mat.cols(), mat.rows(), band);
}
 
// Set a block of data for given x/y position and band from EMatrix2D<T>
template<typename T>
void isce3::io::Raster::setBlock(isce3::core::EMatrix2D<T>& mat, size_t xidx, size_t yidx,
                                    size_t band) {
      setBlock(mat.data(), xidx, yidx, mat.cols(), mat.rows(), band);
}

inline void isce3::io::Raster::setGeoTransform(double* arr)
{
    int status = _dataset->SetGeoTransform(arr);
    if (status != 0) {
        throw isce3::except::RuntimeError(ISCE_SRCINFO(), "Could not set GDAL GeoTransform");
    }
}

inline void isce3::io::Raster::setGeoTransform(std::vector<double>& arr)
{
    if (arr.size() != 6) {
        throw isce3::except::LengthError(ISCE_SRCINFO(), "vector is not of size 6");
    }
    setGeoTransform(arr.data());
}

inline void isce3::io::Raster::setGeoTransform(std::valarray<double>& arr)
{
    if (arr.size() != 6) {
        throw isce3::except::LengthError(ISCE_SRCINFO(), "valarray is not of size 6");
    }
    setGeoTransform(&arr[0]);
}

inline void isce3::io::Raster::getGeoTransform(double *arr) const
{
    int status = _dataset->GetGeoTransform(arr);
    if (status != 0) {
        throw isce3::except::RuntimeError(ISCE_SRCINFO(), "Could not fetch GDAL GeoTransform");
    }
}

inline void isce3::io::Raster::getGeoTransform(std::vector<double>& arr) const
{
    if (arr.size() != 6) {
        throw isce3::except::LengthError(ISCE_SRCINFO(), "vector is not of size 6");
    }
    getGeoTransform(arr.data());
}

inline void isce3::io::Raster::getGeoTransform(std::valarray<double>& arr) const
{
    if (arr.size() != 6) {
        throw isce3::except::LengthError(ISCE_SRCINFO(), "valarray is not of size 6");
    }
    getGeoTransform(&arr[0]);
}
 
//Return western edge
inline double isce3::io::Raster::x0() const
{
    double trans[6];
    getGeoTransform(trans);
    return trans[0];
}
 
//Return northern edge
inline double isce3::io::Raster::y0() const
{
    double trans[6];
    getGeoTransform(trans);
    return trans[3];
}
 
//Return EW pixel spacing
inline double isce3::io::Raster::dx() const
{
    double trans[6];
    getGeoTransform(trans);
    return trans[1];
}
 
//Return NS pixel spacing
inline double isce3::io::Raster::dy() const
{
    double trans[6];
    getGeoTransform(trans);
    return trans[5];
}