Skip to content
Snippets Groups Projects
Commit 28618278 authored by Tom Fischer's avatar Tom Fischer
Browse files

improved class GeoLib::Raster

parent f59b8619
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
GeoLib/Raster.cpp
+ 144
− 15
View file @ 28618278
......@@ -10,13 +10,47 @@
* Created on 2011-09-07 by Thomas Fischer
*/
#include <fstream>
#include "Raster.h"
// BaseLib
#include "StringTools.h"
namespace GeoLib {
Raster::Raster(double cell_size, double no_data_val) :
_cell_size(cell_size), _no_data_val(no_data_val)
Raster::Raster(std::size_t n_cols, std::size_t n_rows, double xllcorner, double yllcorner,
double cell_size, double no_data_val, double* raster_data) :
_n_cols(n_cols), _n_rows(n_rows), _ll_pnt(xllcorner, yllcorner, 0.0),
_cell_size(cell_size), _no_data_val(no_data_val), _raster_data(raster_data)
{}
void Raster::refineRaster(std::size_t n_cols, std::size_t n_rows)
{
if (n_rows <= _n_rows || n_cols <= _n_cols) return;
std::size_t row_blk_size(n_rows / _n_rows);
std::size_t col_blk_size(n_cols / _n_cols);
double *new_raster_data(new double[n_rows*n_cols]);
for (std::size_t row(0); row<_n_rows; row++) {
for (std::size_t col(0); col<_n_cols; col++) {
for (std::size_t new_row(row*row_blk_size); new_row<(row+1)*row_blk_size; new_row++) {
for (std::size_t new_col(col*row_blk_size); new_col<(col+1)*col_blk_size; new_col++) {
new_raster_data[new_row*n_cols+new_col] = _raster_data[row*_n_cols+col];
}
}
}
}
std::swap(_raster_data, new_raster_data);
delete [] new_raster_data;
}
Raster::~Raster()
{
if (_raster_data != NULL)
delete [] _raster_data;
}
void Raster::setCellSize(double cell_size)
......@@ -29,43 +63,138 @@ void Raster::setNoDataVal (double no_data_val)
_no_data_val = no_data_val;
}
double* Raster::getRasterFromSurface(Surface const& sfc, size_t &n_x_pnts, size_t &n_y_pnts) const
void Raster::writeRasterAsASC(std::ostream &os) const
{
// write header
os << "ncols " << _n_rows << std::endl;
os << "nrows " << _n_cols << std::endl;
os << "xllcorner " << _ll_pnt[0] << std::endl;
os << "yllcorner " << _ll_pnt[1] << std::endl;
os << "cellsize " << _cell_size << std::endl;
os << "NODATA_value " << _no_data_val << std::endl;
// write data
for (unsigned row(0); row<_n_rows; row++) {
for (unsigned col(0); col<_n_cols; col++) {
os << _raster_data[(_n_rows-row-1)*_n_cols+col] << " ";
}
os << std::endl;
}
}
Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, double no_data_val)
{
Point const& ll (sfc.getAABB().getMinPoint());
Point const& ur (sfc.getAABB().getMaxPoint());
n_x_pnts = static_cast<size_t>(fabs(ur[0]-ll[0]) / _cell_size)+1;
n_y_pnts = static_cast<size_t>(fabs(ur[1]-ll[1]) / _cell_size)+1;
std::size_t n_cols = static_cast<size_t>(fabs(ur[0]-ll[0]) / cell_size)+1;
std::size_t n_rows = static_cast<size_t>(fabs(ur[1]-ll[1]) / cell_size)+1;
const size_t n_triangles (sfc.getNTriangles());
double *z_vals (new double[n_x_pnts*n_y_pnts]);
double *z_vals (new double[n_cols*n_rows]);
if (!z_vals) {
std::cout << "DEBUG: CreateRaster::getRaster " << n_x_pnts << " x " << n_y_pnts << " to big" << std::endl;
std::cout << "DEBUG: CreateRaster::getRaster " << n_cols << " x " << n_rows << " to big" << std::endl;
return NULL;
}
size_t k(0);
for (size_t r(0); r < n_x_pnts; r++) {
for (size_t c(0); c < n_y_pnts; c++) {
const double test_pnt[3] = { ll[0] + r*_cell_size, ll[1] + c*_cell_size, 0};
for (size_t r(0); r < n_cols; r++) {
for (size_t c(0); c < n_rows; c++) {
const double test_pnt[3] = { ll[0] + r*cell_size, ll[1] + c*cell_size, 0};
for (k=0; k<n_triangles; k++) {
if (sfc[k]->containsPoint2D(test_pnt)) {
Triangle const * const tri (sfc[k]);
// compute coefficients c0, c1, c2 for the plane f(x,y) = c0 x + c1 y + c2
double coeff[3] = {0.0, 0.0, 0.0};
GeoLib::getPlaneCoefficients(*tri, coeff);
z_vals[r*n_y_pnts+c] = coeff[0] * test_pnt[0] + coeff[1] * test_pnt[1] + coeff[2];
z_vals[r*n_rows+c] = coeff[0] * test_pnt[0] + coeff[1] * test_pnt[1] + coeff[2];
break;
}
}
if (k==n_triangles) {
z_vals[r*n_y_pnts+c] = _no_data_val;
z_vals[r*n_rows+c] = no_data_val;
}
}
}
return z_vals;
return new Raster(n_cols, n_rows, ll[0], ll[1], cell_size, -9999, z_vals);
}
Raster::~Raster()
{}
Raster* Raster::getRasterFromASCFile(std::string const& fname)
{
std::ifstream in(fname.c_str());
if (!in.is_open()) {
std::cout << "Raster::getRasterFromASCFile() - Could not open file..." << fname << std::endl;
return NULL;
}
// header information
std::size_t n_cols(0), n_rows(0);
double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), no_data_val(-9999);
if (readASCHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, no_data_val)) {
double* values = new double[n_cols*n_rows];
std::string s;
// read the data into the double-array
for (size_t j(0); j < n_rows; ++j) {
size_t idx ((n_rows - j - 1) * n_cols);
for (size_t i(0); i < n_cols; ++i) {
in >> s;
values[idx+i] = strtod(BaseLib::replaceString(",", ".", s).c_str(),0);
}
}
in.close();
return new Raster(n_cols, n_rows, xllcorner, yllcorner,
cell_size, no_data_val, values);
} else {
std::cout << "Raster::getRasterFromASCFile() - could not read header of file " << fname << std::endl;
return NULL;
}
}
bool Raster::readASCHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val)
{
std::string tag, value;
in >> tag;
if (tag.compare("ncols") == 0) {
in >> value;
n_cols = atoi(value.c_str());
} else return false;
in >> tag;
if (tag.compare("nrows") == 0) {
in >> value;
n_rows = atoi(value.c_str());
} else return false;
in >> tag;
if (tag.compare("xllcorner") == 0) {
in >> value;
xllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
in >> tag;
if (tag.compare("yllcorner") == 0) {
in >> value;
yllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
in >> tag;
if (tag.compare("cellsize") == 0) {
in >> value;
cell_size = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
in >> tag;
if (tag.compare("NODATA_value") == 0) {
in >> value;
no_data_val = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
return true;
}
} // end namespace GeoLib
GeoLib/Raster.h
+ 21
− 4
View file @ 28618278
......@@ -19,14 +19,31 @@ namespace GeoLib {
class Raster {
public:
Raster(double cell_size=1, double no_data_val=9999);
void setCellSize(double cell_size);
void setNoDataVal (double no_data_val);
double* getRasterFromSurface (Surface const& sfc, std::size_t &n_x_pnts, std::size_t &n_y_pnts) const;
Raster(std::size_t n_cols, std::size_t n_rows, double xllcorner, double yllcorner,
double cell_size = 1, double no_data_val = -9999, double* raster_data = NULL);
std::size_t getNCols() const { return _n_cols; }
std::size_t getNRows() const { return _n_rows; }
void refineRaster(std::size_t n_cols, std::size_t n_rows);
double const* getRasterData() const { return _raster_data; }
virtual ~Raster();
void writeRasterAsASC(std::ostream &os) const;
static Raster* getRasterFromASCFile(std::string const& fname);
static Raster* getRasterFromSurface(Surface const& sfc, double cell_size, double no_data_val = -9999);
private:
static bool readASCHeader(std::ifstream &in, size_t &n_cols, std::size_t &n_rows,
double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val);
void setCellSize(double cell_size);
void setNoDataVal (double no_data_val);
std::size_t _n_cols;
std::size_t _n_rows;
GeoLib::Point _ll_pnt;
double _cell_size;
double _no_data_val;
double* _raster_data;
};
}
......
Gui/mainwindow.cpp
+ 15
− 37
View file @ 28618278
......@@ -19,6 +19,9 @@
#include "MshModel.h"
#include "StationTreeModel.h"
// GeoLib
#include "Raster.h"
//dialogs
#include "CondFromRasterDialog.h"
#include "ConditionWriterDialog.h"
......@@ -1130,16 +1133,23 @@ void MainWindow::FEMTestStart()
std::string path("/mnt/visdata/tom/data/Influins/Mapping/TestCase-10x10-100x100/");
// read properties from asc file
std::string fname_asc(path+"Kf-Raster-10x10-100x100.asc");
double x0(0.0), y0(0.0), delta(0.0);
unsigned n_cols(0), n_rows(0);
float* img_data(VtkRaster::loadDataFromASC(fname_asc, x0, y0, n_cols, n_rows, delta));
GeoLib::Raster* raster(GeoLib::Raster::getRasterFromASCFile(fname_asc));
// // write new asc file
// std::ofstream out(path+"Kf-Raster-10x10-100x100.asc");
// raster->refine(raster->getNCols() * 10, raster->getNRows() * 10);
// raster->writeRasterAsASC(out);
// out.close();
double const*const raster_data(raster->getRasterData());
unsigned n_cols(raster->getNCols()), n_rows(raster->getNRows());
std::vector<double> src_properties(n_cols*n_rows);
for (unsigned row(0); row<n_rows; row++) {
for (unsigned col(0); col<n_cols; col++) {
src_properties[row*n_cols+col] = img_data[2*(row*n_cols+col)];
src_properties[row*n_cols+col] = raster_data[2*(row*n_cols+col)];
}
}
delete [] img_data;
delete raster;
{
double src_mean_value(src_properties[0]);
......@@ -1157,38 +1167,6 @@ void MainWindow::FEMTestStart()
std::cout << "variance of source: " << src_varianz << std::endl;
}
// unsigned new_n_rows(100), new_n_cols(100);
// unsigned row_blk_size(new_n_rows / n_rows);
// unsigned col_blk_size(new_n_cols / n_cols);
// std::vector<double> new_src_properties(new_n_cols*new_n_rows);
// for (unsigned row(0); row<n_rows; row++) {
// for (unsigned col(0); col<n_cols; col++) {
// for (unsigned new_row(row*row_blk_size); new_row<(row+1)*row_blk_size; new_row++) {
// for (unsigned new_col(col*row_blk_size); new_col<(col+1)*col_blk_size; new_col++) {
// new_src_properties[new_row*new_n_cols+new_col] = src_properties[row*n_cols+col];
// }
// }
// }
// }
//
// // write new asc file
// std::ofstream out(path+"Kf-Raster-10x10-100x100.asc");
// // write header
// out << "ncols " << new_n_rows << std::endl;
// out << "nrows " << new_n_cols << std::endl;
// out << "xllcorner " << x0 << std::endl;
// out << "yllcorner " << y0 << std::endl;
// out << "cellsize " << delta / row_blk_size << std::endl;
// out << "NODATA_value -9999" << std::endl;
// // write data
// for (unsigned new_row(0); new_row<new_n_rows; new_row++) {
// for (unsigned new_col(0); new_col<new_n_cols; new_col++) {
// out << new_src_properties[(new_n_rows-new_row-1)*new_n_cols+new_col] << " ";
// }
// out << std::endl;
// }
// out.close();
std::vector<size_t> src_perm(n_cols*n_rows);
for (size_t k(0); k<n_cols*n_rows; k++) src_perm[k] = k;
BaseLib::Quicksort<double>(src_properties, 0, n_cols*n_rows, src_perm);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment