Program Listing for File Map.h¶
↰ Return to documentation for file (necsim/Map.h
)
// This file is part of necsim project which is released under MIT license.
// See file **LICENSE.txt** or visit https://opensource.org/licenses/MIT) for full license details
#ifndef MAP_H
#define MAP_H
#ifdef with_gdal
#include <string>
#include <cstring>
#include <sstream>
#include <gdal_priv.h>
#include <cpl_conv.h> // for CPLMalloc()
#include <sstream>
#include "Logging.h"
#include "Matrix.h"
#include "custom_exceptions.h"
#include "cpl_custom_handler.h"
using namespace std;
#ifdef DEBUG
#include "custom_exceptions.h"
#endif // DEBUG
namespace necsim
{
template<class T>
class Map : public virtual Matrix<T>
{
protected:
shared_ptr<GDALDataset*> po_dataset;
shared_ptr<GDALRasterBand*> po_band;
unsigned long block_x_size, block_y_size;
double no_data_value;
string file_name;
GDALDataType gdal_data_type;
// Contains the error object to check for any gdal issues
CPLErr cpl_error;
double upper_left_x, upper_left_y, x_res, y_res;
using Matrix<T>::matrix;
using Matrix<T>::num_cols;
using Matrix<T>::num_rows;
bool cpl_error_set;
public:
using Matrix<T>::setSize;
using Matrix<T>::getCols;
using Matrix<T>::getRows;
using Matrix<T>::get;
using Matrix<T>::begin;
using Matrix<T>::end;
using Matrix<T>::operator*;
using Matrix<T>::operator/;
using Matrix<T>::operator+;
using Matrix<T>::operator-;
using Matrix<T>::operator-=;
using Matrix<T>::operator+=;
using Matrix<T>::operator*=;
using Matrix<T>::operator/=;
using Matrix<T>::operator=;
Map() : Matrix<T>(0, 0), po_dataset(nullptr), po_band(nullptr), block_x_size(0), block_y_size(0),
no_data_value(0.0), file_name(""), gdal_data_type(GDT_Unknown), cpl_error(CE_None), upper_left_x(0.0),
upper_left_y(0.0), x_res(0.0), y_res(0.0), cpl_error_set(false)
{
setCPLErrorHandler();
GDALAllRegister();
removeCPLErrorHandler();
}
~Map() override
{
close();
removeCPLErrorHandler();
}
void setCPLErrorHandler()
{
if(!cpl_error_set)
{
cpl_error_set = true;
CPLPushErrorHandler(cplNecsimCustomErrorHandler);
}
}
void removeCPLErrorHandler()
{
if(cpl_error_set)
{
cpl_error_set = false;
CPLPopErrorHandler();
}
}
void open(const string &filename_in)
{
if(!po_dataset)
{
file_name = filename_in;
setCPLErrorHandler();
// GDALDataset * tmp_ptr = (GDALDataset*)GDALOpen(file_name.c_str(), GA_ReadOnly);
po_dataset = make_shared<GDALDataset*>((GDALDataset*)GDALOpen(file_name.c_str(), GA_ReadOnly));
removeCPLErrorHandler();
}
else
{
throw FatalException("File already open at " + file_name);
}
if(!po_dataset)
{
string s = "File " + file_name + " not found.";
throw runtime_error(s);
}
}
void open()
{
open(file_name);
}
bool isOpen()
{
return po_dataset != nullptr;
}
void close()
{
if(po_dataset)
{
setCPLErrorHandler();
if(po_dataset.use_count() == 1)
{
GDALClose(*po_dataset);
}
removeCPLErrorHandler();
po_dataset = nullptr;
po_band = nullptr;
}
}
void getRasterBand()
{
setCPLErrorHandler();
if((*po_dataset)->GetRasterCount() < 1)
{
stringstream ss;
ss << "No raster band detected in " << file_name << endl;
throw FatalException(ss.str());
}
po_band = make_shared<GDALRasterBand*>((GDALRasterBand*)(*po_dataset)->GetRasterBand(1));
removeCPLErrorHandler();
if(po_band == nullptr)
{
stringstream ss;
ss << "Could not open po_band in " << file_name << ": returned a nullptr" << endl;
throw FatalException(ss.str());
}
}
void getBlockSizes()
{
setCPLErrorHandler();
block_x_size = static_cast<unsigned long>((*po_dataset)->GetRasterXSize());
block_y_size = static_cast<unsigned long>((*po_dataset)->GetRasterYSize());
removeCPLErrorHandler();
}
void getMetaData()
{
#ifdef DEBUG
if(po_dataset == nullptr)
{
throw FatalException("po_dataset is nullptr. This is likely a problem with gdal. Please report this bug.");
}
if(po_band == nullptr)
{
throw FatalException("po_band is nullptr. This is likely a problem with gdal. Please report this bug.");
}
#endif // DEBUG
setCPLErrorHandler();
try
{
int pbSuccess;
no_data_value = (*po_band)->GetNoDataValue(&pbSuccess);
if(!pbSuccess)
{
no_data_value = 0.0;
}
}
catch(out_of_range &out_of_range1)
{
no_data_value = 0.0;
}
stringstream ss;
ss << "No data value is: " << no_data_value << endl;
writeInfo(ss.str());
// Check sizes match
gdal_data_type = (*po_band)->GetRasterDataType();
const static double default_values[6] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
double geoTransform[6] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
writeInfo("Getting geo transform...");
cpl_error = (*po_dataset)->GetGeoTransform(geoTransform);
if(cpl_error >= CE_Warning)
{
ss.str("");
ss << "cpl error detected greater than or equal to warning level." << endl;
writeInfo(ss.str());
string message = "No transform present in dataset for " + file_name;
cplNecsimCustomErrorHandler(cpl_error, 6, message.c_str());
CPLErrorReset();
cpl_error = CE_None;
std::copy(std::begin(default_values), std::end(default_values), std::begin(geoTransform));
}
writeInfo("done.\n");
ss.str("");
ss << "Affine transform is " << geoTransform[0] << ", " << geoTransform[1] << ", " << geoTransform[2]
<< ", ";
ss << geoTransform[3] << ", " << geoTransform[4] << ", " << geoTransform[5] << endl;
writeInfo(ss.str());
upper_left_x = geoTransform[0];
upper_left_y = geoTransform[3];
x_res = geoTransform[1];
y_res = -geoTransform[5];
removeCPLErrorHandler();
// checkTifImportFailure();
#ifdef DEBUG
printMetaData();
#endif // DEBUG
}
#ifdef DEBUG
void printMetaData()
{
setCPLErrorHandler();
stringstream ss;
const char* dt_name = GDALGetDataTypeName(gdal_data_type);
ss << "Filename: " << file_name << endl;
writeLog(10, ss.str());
ss.str("");
ss << "data type: " << gdal_data_type << "(" << dt_name << ")" << endl;
writeLog(10, ss.str());
ss.str("");
ss << "Geo-transform (ulx, uly, x res, y res): " << upper_left_x << ", " << upper_left_y << ", ";
ss << x_res << ", " << y_res << ", " << endl;
writeLog(10, ss.str());
ss.str("");
ss << "No data value: " << no_data_value << endl;
writeLog(10, ss.str());
removeCPLErrorHandler();
}
#endif //DEBUG
double getUpperLeftX() const
{
return upper_left_x;
}
double getUpperLeftY() const
{
return upper_left_y;
}
string getFileName() const
{
return file_name;
}
void import(const string &filename) override
{
if(!importTif(filename))
{
Matrix<T>::import(filename);
}
}
bool importTif(const string &filename)
{
if(filename.find(".tif") != string::npos)
{
setCPLErrorHandler();
stringstream ss;
ss << "Importing " << filename << " " << endl;
writeInfo(ss.str());
open(filename);
getRasterBand();
getBlockSizes();
getMetaData();
// If the sizes are 0 then use the raster sizes
if(num_cols == 0 || num_rows == 0)
{
setSize(block_y_size, block_x_size);
}
// Check sizes
if((num_cols != block_x_size || num_rows != block_y_size) || num_cols == 0 || num_rows == 0)
{
stringstream stringstream1;
stringstream1 << "Raster data size does not match inputted dimensions for " << filename
<< ". Using raster sizes." << endl;
stringstream1 << "Old dimensions: " << num_cols << ", " << num_rows << endl;
stringstream1 << "New dimensions: " << block_x_size << ", " << block_y_size << endl;
writeWarning(stringstream1.str());
setSize(block_y_size, block_x_size);
}
// Check the data types are support
const char* dt_name = GDALGetDataTypeName(gdal_data_type);
if(gdal_data_type == 0 || gdal_data_type > 7)
{
removeCPLErrorHandler();
throw FatalException("Data type of " + string(dt_name) + " is not supported.");
}
#ifdef DEBUG
if(sizeof(T) * 8 != gdal_data_sizes[gdal_data_type])
{
stringstream ss2;
ss2 << "Object data size: " << sizeof(T) * 8 << endl;
ss2 << "Tif data type: " << dt_name << ": " << gdal_data_sizes[gdal_data_type] << " bytes" << endl;
ss2 << "Tif data type does not match object data size in " << file_name << endl;
writeWarning(ss2.str());
}
#endif
// Use the overloaded method for importing between types
internalImport();
writeInfo("done.\n");
removeCPLErrorHandler();
return true;
}
return false;
}
bool openOffsetMap(Map &offset_map)
{
bool opened_here = false;
if(!offset_map.isOpen())
{
opened_here = true;
offset_map.open();
}
offset_map.getRasterBand();
offset_map.getMetaData();
return opened_here;
}
void closeOffsetMap(Map &offset_map, const bool &opened_here)
{
if(opened_here)
{
offset_map.close();
}
}
void calculateOffset(Map &offset_map, long &offset_x, long &offset_y)
{
auto opened_here = openOffsetMap(offset_map);
offset_x = static_cast<long>(round((upper_left_x - offset_map.upper_left_x) / x_res));
offset_y = static_cast<long>(round((offset_map.upper_left_y - upper_left_y) / y_res));
closeOffsetMap(offset_map, opened_here);
}
unsigned long roundedScale(Map &offset_map)
{
auto opened_here = openOffsetMap(offset_map);
closeOffsetMap(offset_map, opened_here);
return static_cast<unsigned long>(floor(offset_map.x_res / x_res));
}
void internalImport()
{
writeWarning(
"\nNo type detected for Map type. Attempting default importing (potentially undefined behaviour).");
defaultImport();
}
void defaultImport()
{
setCPLErrorHandler();
unsigned int number_printed = 0;
for(uint32_t j = 0; j < num_rows; j++)
{
printNumberComplete(j, number_printed);
cpl_error = (*po_band)->RasterIO(GF_Read, 0, j, static_cast<int>(block_x_size), 1, &matrix[j * num_cols],
static_cast<int>(block_x_size), 1, gdal_data_type, 0, 0);
checkTifImportFailure();
// Now convert the no data values to 0
for(uint32_t i = 0; i < num_cols; i++)
{
if(get(j, i) == no_data_value)
{
get(j, i) = 0;
}
}
}
removeCPLErrorHandler();
}
void importFromDoubleAndMakeBool()
{
setCPLErrorHandler();
writeInfo("\nConverting from double to boolean.\n");
unsigned int number_printed = 0;
// create an empty row of type float
double* t1;
t1 = (double*) CPLMalloc(sizeof(double) * num_cols);
// import the data a row at a time, using our template row.
for(uint32_t j = 0; j < num_rows; j++)
{
printNumberComplete(j, number_printed);
cpl_error = (*po_band)->RasterIO(GF_Read, 0, j, static_cast<int>(block_x_size), 1, &t1[0],
static_cast<int>(block_x_size), 1, GDT_Float64, 0, 0);
checkTifImportFailure();
// now copy the data to our Map, converting float to int. Round or floor...? hmm, floor?
for(unsigned long i = 0; i < num_cols; i++)
{
if(t1[i] == no_data_value)
{
this->setValue(j, i, false);
}
else
{
this->setValue(j, i, t1[i] >= 0.5);
}
}
}
CPLFree(t1);
removeCPLErrorHandler();
}
template<typename T2> void importUsingBuffer(GDALDataType dt_buff)
{
setCPLErrorHandler();
stringstream ss;
ss << "\nUsing buffer of type " << GDALGetDataTypeName(dt_buff) << " into array of dimensions ";
ss << num_cols << " by " << num_rows << " with block size " << block_x_size << ", " << block_y_size << endl;
writeInfo(ss.str());
unsigned int number_printed = 0;
// create an empty row of type float
T2* t1;
t1 = (T2*) CPLMalloc(sizeof(T2) * num_cols);
if(CPLGetLastErrorNo() != CE_None)
{
stringstream os;
os << "Error thrown by CPLMalloc: " << CPLGetLastErrorType() << ": " << CPLGetLastErrorMsg() << endl;
throw FatalException(os.str());
}
// import the data a row at a time, using our template row.
auto int_block_x_size = static_cast<int>(block_x_size);
#ifdef DEBUG
if(sizeof(T2) * 8 != GDALGetDataTypeSize(dt_buff))
{
stringstream ss0;
ss0 << "Size of template (" << sizeof(T2) << ") does not equal size of gdal buffer (";
ss0 << GDALGetDataTypeSize(dt_buff) << "). Please report this bug." << endl;
throw FatalException(ss0.str());
}
if(t1 == nullptr)
{
stringstream ss1;
ss1 << "CPL malloc could not acquire " << sizeof(T2) * num_cols << " of space and returned a nullptr.";
ss1 << " Please check that your install of gdal is fully functional, and otherwise report this bug."
<< endl;
throw FatalException(ss1.str());
}
if(static_cast<unsigned long>(int_block_x_size) != block_x_size)
{
stringstream ss2;
ss2 << "Error in integer conversion - please report this bug: " << int_block_x_size << " != ";
ss2 << block_x_size << endl;
throw FatalException(ss2.str());
}
if(po_band == nullptr)
{
CPLFree(t1);
throw FatalException("po_band is nullptr during import using buffer - please report this bug.");
}
#endif // DEBUG
for(uint32_t j = 0; j < num_rows; j++)
{
printNumberComplete(j, number_printed);
cpl_error = (*po_band)->RasterIO(GF_Read, 0, j, int_block_x_size, 1, t1, int_block_x_size, 1, dt_buff, 0,
0);
checkTifImportFailure();
for(unsigned long i = 0; i < num_cols; i++)
{
if(t1[i] == no_data_value)
{
this->setValue(j, i, static_cast<T>(0));
}
else
{
this->setValue(j, i, static_cast<T>(t1[i]));
}
}
}
CPLFree(t1);
removeCPLErrorHandler();
}
void printNumberComplete(const uint32_t &j, unsigned int &number_printed)
{
double dComplete = ((double) j / (double) num_rows) * 20;
if(number_printed < dComplete)
{
stringstream os;
os << "\rImporting " << file_name << " ";
number_printed = 0;
while(number_printed < dComplete)
{
os << ".";
number_printed++;
}
os << flush;
writeInfo(os.str());
}
}
void checkTifImportFailure()
{
if(cpl_error >= CE_Warning && !cpl_error_set)
{
stringstream ss;
ss << "\nCPL error thrown during import of " << file_name << endl;
cplNecsimCustomErrorHandler(cpl_error, 3, ss.str().c_str());
CPLErrorReset();
cpl_error = CE_None;
}
}
friend ostream &operator>>(ostream &os, const Map &m)
{
return Matrix<T>::writeOut(os, m);
}
friend istream &operator<<(istream &is, Map &m)
{
return Matrix<T>::readIn(is, m);
}
Map &operator=(const Map &m)
{
Matrix<T>::operator=(m);
this->po_dataset = m.po_dataset;
this->po_band = m.po_band;
this->block_x_size = m.block_x_size;
this->block_y_size = m.block_y_size;
this->no_data_value = m.no_data_value;
this->file_name = m.file_name;
this->gdal_data_type = m.gdal_data_type;
this->cpl_error = m.cpl_error;
this->upper_left_x = m.upper_left_x;
this->upper_left_y = m.upper_left_y;
this->x_res = m.x_res;
this->y_res = m.y_res;
this->cpl_error_set = m.cpl_error_set;
return *this;
}
};
template<> inline void Map<bool>::internalImport()
{
if(gdal_data_type <= 5)
{
// Then the tif file type is an int/byte
// we can just import as it is
importUsingBuffer<uint8_t>(GDT_Byte);
}
else
{
// Conversion from double to boolean
importFromDoubleAndMakeBool();
}
}
template<> inline void Map<int8_t>::internalImport()
{
importUsingBuffer<int16_t>(GDT_Int16);
}
template<> inline void Map<uint8_t>::internalImport()
{
gdal_data_type = GDT_Byte;
defaultImport();
}
template<> inline void Map<int16_t>::internalImport()
{
gdal_data_type = GDT_Int16;
defaultImport();
}
template<> inline void Map<uint16_t>::internalImport()
{
gdal_data_type = GDT_UInt16;
defaultImport();
}
template<> inline void Map<int32_t>::internalImport()
{
gdal_data_type = GDT_Int32;
defaultImport();
}
template<> inline void Map<uint32_t>::internalImport()
{
gdal_data_type = GDT_UInt32;
defaultImport();
}
template<> inline void Map<float>::internalImport()
{
gdal_data_type = GDT_Float32;
defaultImport();
}
template<> inline void Map<double>::internalImport()
{
gdal_data_type = GDT_Float64;
defaultImport();
}
#endif // with_gdal
}
#endif //MAP_H