Commit a3e53b6a authored by mjg226's avatar mjg226
Browse files

bug fixing and updating asciigrid for GDAL usage

parent e6b8d751
......@@ -260,8 +260,6 @@ CAAPI_add_source_dir(impls\\\\specific "${CAAPI_SPECIFIC_IMPL_DIR}/*.ca"
set(CAAPI_IMPL_INCLUDE_DIRS)
set(CAAPI_IMPL_LIBRARIES)
# Add the implementation dir to the includ dirs
list(APPEND CAAPI_IMPL_INCLUDE_DIRS ${CAAPI_SPECIFIC_IMPL_DIR})
# Execute the module for the specific implementation
if(EXISTS "${CAAPI_SPECIFIC_IMPL_DIR}/extra/CAAPI_impl.cmake")
......@@ -271,6 +269,10 @@ if(EXISTS "${CAAPI_SPECIFIC_IMPL_DIR}/extra/CAAPI_impl.cmake")
CAAPI_impl("${CAAPI_SPECIFIC_IMPL_DIR}" CAAPI_IMPL_INCLUDE_DIRS CAAPI_IMPL_LIBRARIES)
endif()
# Add the implementation dir to the includ dirs
list(APPEND CAAPI_IMPL_INCLUDE_DIRS ${CAAPI_SPECIFIC_IMPL_DIR})
# Execute the module that manage the common ca for the specific implementation
if(EXISTS "${CAAPI_SPECIFIC_IMPL_DIR}/extra/CAAPI_ca.cmake")
......@@ -286,6 +288,8 @@ add_custom_target(caapi_do_impl_headers ALL DEPENDS ${CAAPI_CAHEADERS} ${CAAPI_H
# Add the needed directories into the compilation
include_directories(${CAAPI_IMPL_INCLUDE_DIRS})
message(STATUS "Include dir: ${CAAPI_IMPL_INCLUDE_DIRS}")
########################## CAPI DEVICE common files ##########################
......
......@@ -30,7 +30,7 @@ else()
add_definitions( -DCAAPI_IMPL_VERSION=${CAAPI_IMPL_VERSION} )
endif()
if (WIN32)
if (WIN32 AND NOT MINGW)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /openmp")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /openmp")
else()
......
/*
Copyright (c) 2013 Centre for Water Systems,
Copyright (c) 2016 Centre for Water Systems,
University of Exeter
Permission is hereby granted, free of charge, to any person obtaining a copy
......@@ -23,8 +23,6 @@ THE SOFTWARE.
*/
#ifndef _CA_ESRI_ASCIIGRID_HPP_
#define _CA_ESRI_ASCIIGRID_HPP_
//! \file AsciiGrid.hpp
......@@ -33,944 +31,9 @@ THE SOFTWARE.
//! \author Michele Guidolin, University of Exeter,
//! \author Mike Gibson, University of Exeter,
//! contact: m.j.gibson [at] exeter.ac.uk
//! \date 2012-01
//! \date 2016-01
#include "AsciiGrid.hpp"
#include "../simple/ESRI_ASCIIGrid.hpp"
namespace CA {
template<typename T>
class ESRI_ASCIIGrid
: public CA::AsciiGridGeneral<T>
{
public:
//public ESRI_ASCIIGrid<T>() : CA::AsciiGridGeneral<T>(){
//};
//! Read token and value pair from a line of the header of an ARC/INFO
//! ASCII GRID file.
template<typename A>
inline
void readAsciiGridHeaderLine(const std::string& filename, std::ifstream& infile, A& value,
const char* check, bool substring = false)
{
std::string token;
std::string strvalue;
// Read the token and value. Check that the token is
// the desire one. Add the value into grid.ncols.
if (!(infile >> token) || !(infile >> strvalue))
throw std::runtime_error(std::string("Error reading the file: ") + filename);
if (!CA::compareCaseInsensitive(token, check, substring))
throw std::runtime_error(std::string("Error not an ARC/INFO ASCII GRID file: ") + filename);
if (!CA::fromString(value, strvalue))
throw std::runtime_error(std::string("Error converting the string ") + strvalue
+ std::string(" into a value from the file: ") + filename);
}
//! Read a ARC/INFO ASCII GRID format file
//template<typename T>
inline void readAsciiGrid(const std::string& filename, bool print = false)
{
// Get starting time.
//Clock total_timer;
//std::cout << "READING DATA!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
// Open file
std::ifstream infile(filename.c_str());
// Check if the file was open.
if (!infile)
throw std::runtime_error(std::string("Error opening file: ") + filename);
// Read the header
readAsciiGridHeaderLine<CA::Unsigned>(filename, infile, AsciiGridGeneral<T>::ncols, "ncols", false);
readAsciiGridHeaderLine<CA::Unsigned>(filename, infile, AsciiGridGeneral<T>::nrows, "nrows", false);
readAsciiGridHeaderLine<double>(filename, infile, AsciiGridGeneral<T>::xllcorner, "xllcorner", false);
readAsciiGridHeaderLine<double>(filename, infile, AsciiGridGeneral<T>::yllcorner, "yllcorner", false);
readAsciiGridHeaderLine<double>(filename, infile, AsciiGridGeneral<T>::cellsize, "cellsize", false);
readAsciiGridHeaderLine<T>(filename, infile, AsciiGridGeneral<T>::nodata, "nodata_value", true);
if (print)
{
std::cout << "ncols \t\t" << AsciiGridGeneral<T>::ncols << std::endl;
std::cout << "nrows \t\t" << AsciiGridGeneral<T>::nrows << std::endl;
std::cout << "xllcorner \t" << AsciiGridGeneral<T>::xllcorner << std::endl;
std::cout << "yllcorner \t" << AsciiGridGeneral<T>::yllcorner << std::endl;
std::cout << "cellsize \t" << AsciiGridGeneral<T>::cellsize << std::endl;
std::cout << "nodata \t\t" << AsciiGridGeneral<T>::nodata << std::endl;
std::cout << std::endl;
}
// Allocate data and set the default value to nodata.
AsciiGridGeneral<T>::data.clear();
AsciiGridGeneral<T>::data.resize(AsciiGridGeneral<T>::ncols * AsciiGridGeneral<T>::nrows, AsciiGridGeneral<T>::nodata);
// Read data.
T value;
size_t i = 0;
while (infile >> value && i < AsciiGridGeneral<T>::data.size())
{
if (infile.fail())
throw std::runtime_error(std::string("Error converting a data string into a value"));
AsciiGridGeneral<T>::data[i] = value;
i++;
}
// Close the file.
infile.close();
//std::cout<<"Read run time taken (s) = " << total_timer.millisecond()/1000.0 << std::endl;
}
//! Read a ARC/INFO ASCII GRID format file without loading teh data.
//template<typename T>
inline void readAsciiGridHeader(const std::string& filename, bool print = false)
{
// Open file
std::ifstream infile(filename.c_str());
// Check if the file was open.
if (!infile)
throw std::runtime_error(std::string("Error opening file: ") + filename);
// Read the header
readAsciiGridHeaderLine<CA::Unsigned>(filename, infile, AsciiGridGeneral<T>::ncols, "ncols", false);
readAsciiGridHeaderLine<CA::Unsigned>(filename, infile, AsciiGridGeneral<T>::nrows, "nrows", false);
readAsciiGridHeaderLine<double>(filename, infile, AsciiGridGeneral<T>::xllcorner, "xllcorner", false);
readAsciiGridHeaderLine<double>(filename, infile, AsciiGridGeneral<T>::yllcorner, "yllcorner", false);
readAsciiGridHeaderLine<double>(filename, infile, AsciiGridGeneral<T>::cellsize, "cellsize", false);
readAsciiGridHeaderLine<T>(filename, infile, AsciiGridGeneral<T>::nodata, "nodata_value", true);
if (print)
{
std::cout << "ncols \t\t" << AsciiGridGeneral<T>::ncols << std::endl;
std::cout << "nrows \t\t" << AsciiGridGeneral<T>::nrows << std::endl;
std::cout << "xllcorner \t" << AsciiGridGeneral<T>::xllcorner << std::endl;
std::cout << "yllcorner \t" << AsciiGridGeneral<T>::yllcorner << std::endl;
std::cout << "cellsize \t" << AsciiGridGeneral<T>::cellsize << std::endl;
std::cout << "nodata \t\t" << AsciiGridGeneral<T>::nodata << std::endl;
std::cout << std::endl;
}
// Empty the data
AsciiGridGeneral<T>::data.clear();
// Close the file.
infile.close();
}
//! Write a ARC/INFO ASCII GRID format file
//template<typename T>
inline void writeAsciiGrid(const std::string& filename,
int decimal_places = 6, bool print = false)
{
// Get starting time.
//Clock total_timer;
// Add file extension
String _filename = filename + ".asc";
// Open file
std::ofstream onfile(_filename.c_str());
// Check if the file was open.
if (!onfile)
throw std::runtime_error(std::string("Error opening file: " + _filename));
// Set manipulators
onfile.unsetf(std::ios::floatfield);
onfile.precision(12);
// Write the header
onfile << "ncols \t\t" << AsciiGridGeneral<T>::ncols << std::endl;
onfile << "nrows \t\t" << AsciiGridGeneral<T>::nrows << std::endl;
onfile << "xllcorner \t" << AsciiGridGeneral<T>::xllcorner << std::endl;
onfile << "yllcorner \t" << AsciiGridGeneral<T>::yllcorner << std::endl;
onfile << "cellsize \t" << AsciiGridGeneral<T>::cellsize << std::endl;
onfile << "NODATA_value \t\t" << AsciiGridGeneral<T>::nodata << std::endl;
onfile << std::endl;
// Make sure the decimal places is not zeros.
if (decimal_places <= 0)
decimal_places = 6;
onfile.unsetf(std::ios::floatfield);
onfile << std::fixed;
onfile.precision(decimal_places);
// Loop through the data.
for (unsigned int i = 0; i < AsciiGridGeneral<T>::ncols * AsciiGridGeneral<T>::nrows; ++i)
{
if (i%AsciiGridGeneral<T>::ncols == 0)
onfile << "\n";
onfile << AsciiGridGeneral<T>::data[i] << " ";
}
// End of line.
onfile << std::endl;
// Close the file.
onfile.close();
//std::cout<<"Write run time taken (s) = " << total_timer.millisecond()/1000.0 << std::endl;
}
//------------------------------------------------------------------------------------------------------------------------------
//! Add another ascii gri dof data to this grid of data
inline void addAsciiGrid(AsciiGrid<T>& grid2, bool print = false) {
// find how many cells between bottom left of each terrain.
int colsDiff = 0;
int rowsDiff = 0;
if (print)
std::cout << "BLAH" << std::endl;
// used to be a utility for converting real to unsigned, bewarned linux ppl
colsDiff = ((CA::Unsigned)grid2.xllcorner) - ((CA::Unsigned)AsciiGridGeneral<T>::xllcorner);
rowsDiff = ((CA::Unsigned)grid2.yllcorner) - ((CA::Unsigned)AsciiGridGeneral<T>::yllcorner);
// check for case without alignment?
if (print){
std::cout << "This grid:" << std::endl;
std::cout << "This grid: xllCorner: " << AsciiGridGeneral<T>::xllcorner << std::endl;
std::cout << "This grid: yllCorner: " << AsciiGridGeneral<T>::yllcorner << std::endl;
std::cout << "This grid: nrows : " << AsciiGridGeneral<T>::nrows << std::endl;
std::cout << "This grid: ncols : " << AsciiGridGeneral<T>::ncols << std::endl;
std::cout << "That grid:" << std::endl;
std::cout << "That grid: xllCorner: " << grid2.xllcorner << std::endl;
std::cout << "That grid: yllCorner: " << grid2.yllcorner << std::endl;
std::cout << "That grid: nrows : " << grid2.nrows << std::endl;
std::cout << "That grid: ncols : " << grid2.ncols << std::endl;
std::cout << "Cols Diff : " << colsDiff << std::endl;
std::cout << "Rows Diff : " << rowsDiff << std::endl;
}
// check for overlap in x direction
CA::Unsigned thisXLeft = (CA::Unsigned)AsciiGridGeneral<T>::xllcorner;
CA::Unsigned thisXRight = ((CA::Unsigned)AsciiGridGeneral<T>::xllcorner) + AsciiGridGeneral<T>::ncols;
CA::Unsigned thatXLeft = (CA::Unsigned)grid2.xllcorner;
CA::Unsigned thatXRight = ((CA::Unsigned)grid2.xllcorner) + grid2.ncols;
CA::Unsigned thisYBottom = (CA::Unsigned)AsciiGridGeneral<T>::yllcorner;
CA::Unsigned thisYTop = ((CA::Unsigned)AsciiGridGeneral<T>::yllcorner) + AsciiGridGeneral<T>::nrows;
CA::Unsigned thatYBottom = (CA::Unsigned)grid2.yllcorner;
CA::Unsigned thatYTop = ((CA::Unsigned)grid2.yllcorner) + grid2.nrows;
if (print){
std::cout << "thisXLeft : " << thisXLeft << std::endl;
std::cout << "thatXLeft : " << thatXLeft << std::endl;
std::cout << std::endl;
std::cout << "thisXRight : " << thisXRight << std::endl;
std::cout << "thatXRight : " << thatXRight << std::endl;
std::cout << std::endl;
std::cout << "thisYBottom : " << thisYBottom << std::endl;
std::cout << "thatYBottom : " << thatYBottom << std::endl;
std::cout << std::endl;
std::cout << "thisYTop : " << thisYTop << std::endl;
std::cout << "thatYTop : " << thatYTop << std::endl;
std::cout << std::endl;
}
// if( ( (thatXleft >= thisXleft) && (thatXLeft < thisXRight) ) ||
// ( (thatXRight >= thisXleft) && (thatXRight < thisXRight))
// )
// check for not overlap
if (
((thisXRight < thatXLeft) || (thatXRight < thisXLeft)) // not overlap in x direction...
|| // if no overlap in either or both then not overlapping area
((thisYTop < thatYBottom) || (thatYTop < thisYBottom)) // not overlap in y direction...
)
{
if (print){
std::cout << "No Overlapiing area, so returning from adding therrain data." << std::endl;
}
// no overlapping area, bail out, return message somehow...
return; // bail as nothing to do here
}
// otherwise there is overlap
// calculate the range, i.e. distance of overlaps
CA::Unsigned colsOverlap = 0;
CA::Unsigned rowsOverlap = 0;
if (thatXLeft < thisXLeft) {
// then we start on this grid left most cells.
// next we check to see if that grid totally covers the x direction of this grid
if (thatXRight > thisXRight) {
colsOverlap = AsciiGridGeneral<T>::ncols;
}
else {
// otherwise that grid, end somewhere before this grids x extent.
colsOverlap = thatXRight - thisXLeft;
}
}
else {
//otherwise, that grids left most, starts inside this grids extent
if (thatXRight < thisXRight)
{
// therefore that grids extent is totaly within this grid
colsOverlap = grid2.ncols;
}
else {
// otherwise this grid ends before that one
colsOverlap = thisXRight - thatXLeft;
}
}
if (thatYBottom < thisYBottom) {
// then we start on this grid Bottom most cells.
// next we check to see if that grid totally covers the y direction of this grid
if (thatYTop > thisYTop) {
rowsOverlap = AsciiGridGeneral<T>::nrows;
}
else {
// otherwise that grid, end somewhere before this grids y extent.
rowsOverlap = thatYTop - thisYBottom;
}
}
else {
//otherwise, that grids bottom most, starts inside this grids extent
if (thatYTop < thisYTop)
{
// therefore that grids extent is totaly within this grid
rowsOverlap = grid2.nrows;
}
else {
// otherwise this grid ends before that one
rowsOverlap = thisYTop - thatYBottom;
}
}
// so we can now loop over the area (columns and rows) of the overlapping area.
// NOTE: that geo co-ordinates start at bottom left, and overlap calculations are based on these, therefore
// a conversion to the local co-ordinates which are top left are required.
CA::Unsigned thisXInset = 0;
CA::Unsigned thisYInset = 0;
CA::Unsigned thatXInset = 0;
CA::Unsigned thatYInset = 0;
if (colsDiff >= 0) {
thisXInset = colsDiff;
}
else {
thatXInset = -colsDiff;
}
if (rowsDiff >= 0) {
thisYInset = rowsDiff;
}
else {
thatYInset = -rowsDiff;
}
// now we know the inset due to over lap, based on bottom left.
for (CA::Unsigned col = 0; col < colsOverlap; ++col) {
for (CA::Unsigned row = 0; row < rowsOverlap; ++row) {
// convert to local cordinates (from bottom left to top left)
CA::Unsigned thisGridXCord = col + thisXInset;
CA::Unsigned thisGridYCord = (AsciiGridGeneral<T>::nrows - 1) - (row + thisYInset);
CA::Unsigned thatGridXCord = col + thatXInset;
CA::Unsigned thatGridYCord = (grid2.nrows - 1) - (row + thatYInset);
//T thisDataPoint = AsciiGridGeneral<T>::data[thisGridXCord + (thisGridYCord * AsciiGridGeneral<T>::ncols)];
//T thatDataPoint = grid2.data[thatGridXCord + (thatGridYCord * grid2.ncols)];
//------------------------------------
AsciiGridGeneral<T>::data[thisGridXCord + (thisGridYCord * AsciiGridGeneral<T>::ncols)] =
grid2.data[thatGridXCord + (thatGridYCord * grid2.ncols)];
}// end of looping overlap grids rows
}// end of looping overlap grids coloumns
}// end of addAsciiGrid
//#######################################################################################################################################
//########################################### comparitor functions ########################################################
//#######################################################################################################################################
//--------------------------------------------------------------------------------------------------------------
//! Compare the RMSE of two ascii grids, and return the errors
//! if the same cell from either grid is no data, then no comparison is made
//! note this object is considered the observation (truth) to compare to the other simulation.
inline Errors compareAsciiGrid( AsciiGrid<T>& grid2,T tolerance = 0.0, bool print = false)
{
//
//double result = 0.0;
//int count = 0;
// note needs checking for cell size usuage?!.............
// erros structure used, to reduce multiple loops
Errors result;
result.RMSE = 0.0;
result.RMSE_wet_both = 0.0;
result.RMSE_wet_either = 0.0;
result.meanError = 0.0;
result.meanError_wet_both = 0.0;
result.meanError_wet_either = 0.0;
result.accuracy = 0.0;
result.sensitivity = 0.0;
result.precision = 0.0;
result.precentage = 0.0;
result.nashSutcliffe = 0.0;
result.truePositive = 0;
result.falsePositive = 0;
result.trueNegative = 0;
result.falseNegative = 0;
result.thisWettedCount = 0;
result.thatWettedCount = 0;
result.combinedWettedCount = 0;
result.allDataCells = 0;
result.eitherWetted = 0;
// for nash sutcliffe we need the mean of the observed (truth), i.e. this ascii grids data
double meanObserved = 0.0;
double nashSumSquaredErrors = 0.0;
double nashSumDataVariance = 0.0;
// do a check for size mismatch here.
// due to geo-cor-ordinates, we must check only the region that overlaps
// find how many cells between bottom left of each terrain.
int colsDiff = 0;
int rowsDiff = 0;
// used to be a utility for converting real to unsigned, bewarned linux ppl
colsDiff = ((CA::Unsigned)grid2.xllcorner) - ((CA::Unsigned)AsciiGridGeneral<T>::xllcorner);
rowsDiff = ((CA::Unsigned)grid2.yllcorner) - ((CA::Unsigned)AsciiGridGeneral<T>::yllcorner);
// check for case without alignment?
// check for overlap in x direction
CA::Unsigned thisXLeft = (CA::Unsigned)AsciiGridGeneral<T>::xllcorner;
CA::Unsigned thisXRight = ((CA::Unsigned)AsciiGridGeneral<T>::xllcorner) + AsciiGridGeneral<T>::ncols;
CA::Unsigned thatXLeft = (CA::Unsigned)grid2.xllcorner;
CA::Unsigned thatXRight = ((CA::Unsigned)grid2.xllcorner) + grid2.ncols;
CA::Unsigned thisYBottom = (CA::Unsigned)AsciiGridGeneral<T>::yllcorner;
CA::Unsigned thisYTop = ((CA::Unsigned)AsciiGridGeneral<T>::yllcorner) + AsciiGridGeneral<T>::nrows;
CA::Unsigned thatYBottom = (CA::Unsigned)grid2.yllcorner;
CA::Unsigned thatYTop = ((CA::Unsigned)grid2.yllcorner) + grid2.nrows;
// if( ( (thatXleft >= thisXleft) && (thatXLeft < thisXRight) ) ||
// ( (thatXRight >= thisXleft) && (thatXRight < thisXRight))
// )
// check for not overlap
if(
( (thisXRight < thatXLeft) || (thatXRight < thisXLeft) ) // not overlap in x direction...
|| // if no overlap in either or both then not overlapping area
((thisYTop < thatYBottom) || (thatYTop < thisYBottom)) // not overlap in y direction...
)
{
// no overlapping area, bail out, return message somehow...
}
// otherwise there is overlap