Raster.h

//-*- C++ -*-
//-*- coding: utf-8 -*-
//
// Author: Marco Lavalle
// Original code: Joshua Cohen
// Copyright 2018

#pragma once

#include "forward.h"

#include <complex>
#include <cstdint>
#include <string>
#include <typeindex>
#include <unordered_map>
#include <vector>
#include <valarray>
#include <gdal_priv.h>
#include <gdal_vrt.h>
#include <ogr_spatialref.h>
#include "Constants.h"
#include <isce3/core/Matrix.h>

#include <isce3/io/gdal/Raster.h>

//#include <pyre/journal.h>

class isce3::io::Raster {

  public:

      Raster() {_dataset = nullptr;};

      Raster(const std::string& fname, GDALAccess access = GA_ReadOnly);

      Raster(const std::string& fname, size_t width, size_t length,
             size_t numBands, GDALDataType dtype = isce3::io::defaultGDALDataType,
             const std::string& driverName = isce3::io::defaultGDALDriver);

      Raster(const std::string& fname, size_t width, size_t length, GDALDataType dtype = isce3::io::defaultGDALDataType);

      template<typename T> Raster(isce3::core::Matrix<T> &matrix);

      // Constructor for a 1 band dataset from isce3::core::Matrix<T>::view_type
      template<typename T> Raster(pyre::grid::View<T> &view);

      Raster(const std::string& fname, const Raster& rast);
      template<typename T> Raster(const std::string& fname, const std::vector<T>& buffer,   size_t length) :
            Raster(fname, buffer.size(), length, 1, asGDT<T>) {}
      template<typename T> Raster(const std::string& fname, const std::valarray<T>& buffer, size_t length) :
            Raster(fname, buffer.size(), length, 1, asGDT<T>) {}

      Raster(const std::string& fname, const std::vector<Raster>& rastVec);

      Raster(isce3::io::gdal::Raster&);

      Raster(const Raster&);

      Raster(GDALDataset *inputDataset, bool owner=true);

      inline void initFromPointer(void* ptr, GDALDataType dtype, size_t width, size_t length, size_t pixeloffset, size_t lineoffset);

      ~Raster();

      inline Raster&      operator=(const Raster&);

      inline size_t       length()   const { return _dataset->GetRasterYSize(); }

      inline size_t       width()    const { return _dataset->GetRasterXSize(); }

      inline size_t       numBands() const { return _dataset->GetRasterCount(); }

      inline GDALAccess   access()   const { return _dataset->GetAccess(); }

      inline GDALDataset* dataset()  const { return _dataset; }

      inline void         dataset(GDALDataset* ds) { _dataset=ds; }

      inline GDALDataType dtype(const size_t band=1) const { return _dataset->GetRasterBand(band)->GetRasterDataType(); }

      inline bool         match(const Raster & rast) const { return width()==rast.width() && length()==rast.length(); }

      inline void         open(const std::string &fname, GDALAccess access);

      inline void         addRasterToVRT(const isce3::io::Raster& rast);

      inline void         addBandToVRT(GDALRasterBand *inBand);

      inline void         addRawBandToVRT(const std::string &fname, GDALDataType dtype);
      //void close() { GDALClose( _dataset ); }  // todo: fix segfault conflict with destructor

      //Pixel read/write with buffer passed by reference, optional band index
      template<typename T> void getSetValue(T& buffer, size_t xidz, size_t yidx, size_t band, GDALRWFlag);
      template<typename T> void getValue(T& buffer, size_t xidx, size_t yidx, size_t band = 1);
      template<typename T> void setValue(T& buffer, size_t xidx, size_t yidx, size_t band = 1);

      // Line read/write with raw pointer and width or STL container, optional band index
      template<typename T> void getSetLine(T* buffer, size_t yidx, size_t iowidth, size_t band, GDALRWFlag iodir);
      template<typename T> void getLine(T* buffer,             size_t yidx, size_t iowidth, size_t band = 1);
      template<typename T> void getLine(std::vector<T>& vec,   size_t yidx, size_t band = 1);
      template<typename T> void getLine(std::valarray<T>& arr, size_t yidx, size_t band = 1);
      template<typename T> void setLine(T* buffer,             size_t yidx, size_t iowidth, size_t band = 1);
      template<typename T> void setLine(std::vector<T>& vec,   size_t yidx, size_t band = 1);
      template<typename T> void setLine(std::valarray<T>& arr, size_t yidx, size_t band = 1);

      // 2D block read/write, generic container w/ width or STL container, optional band index
      template<typename T> void getSetBlock(T* buffer,          size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band, GDALRWFlag iodir);
      template<typename T> void getBlock(T* buffer,             size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band = 1);
      template<typename T> void getBlock(std::vector<T>& vec,   size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band = 1);
      template<typename T> void getBlock(std::valarray<T>& arr, size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band = 1);
      template<typename T> void setBlock(T* buffer,             size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band = 1);
      template<typename T> void setBlock(std::vector<T>& vec,   size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band = 1);
      template<typename T> void setBlock(std::valarray<T>& arr, size_t xidx, size_t yidx, size_t iowidth, size_t iolength, size_t band = 1);

      //2D block read/write for Matrix<T>, optional band index
      template<typename T> void getBlock(isce3::core::Matrix<T>& mat, size_t xidx, size_t yidx, size_t band = 1);
      template<typename T> void setBlock(isce3::core::Matrix<T>& mat, size_t xidx, size_t yidx, size_t band = 1);
      //2D block read/write for Matrix<T>, optional band index
      template<typename T> void getSetBlock(pyre::grid::View<T>& view, size_t xidx, size_t yidx, size_t band, GDALRWFlag iodir);
      template<typename T> void    getBlock(pyre::grid::View<T>& view, size_t xidx, size_t yidx, size_t band = 1);
      template<typename T> void    setBlock(pyre::grid::View<T>& view, size_t xidx, size_t yidx, size_t band = 1);

      //Functions to deal with projections and geotransform information
      int getEPSG();
      int setEPSG(int code);
      inline void setGeoTransform(double *arr);
      inline void setGeoTransform(std::vector<double>&);
      inline void setGeoTransform(std::valarray<double>&);
      inline void getGeoTransform(double *) const;
      inline void getGeoTransform(std::vector<double>&) const;
      inline void getGeoTransform(std::valarray<double>&) const;
      //Read only functions for specific elements of GeoTransform
      inline double x0() const;
      inline double y0() const;
      inline double dx() const;
      inline double dy() const;

private:
    GDALDataset * _dataset;
    bool _owner = true;
};

#define ISCE_IO_RASTER_ICC
#include "Raster.icc"
#undef ISCE_IO_RASTER_ICC