<commit_before>/**
Authors:
<NAME> <<EMAIL>>

Maintainers:
Currently maintained by the authors.

This file is part of the util library used by models created at the Institute of
Landscape Systems Analysis at the ZALF.
Copyright (C) 2007-2013, Leibniz Centre for Agricultural Landscape Research (ZALF)

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <fstream>
#include <sstream>
#include <cstddef>
#include <cstdio>
#include <sys/types.h>
#include <sys/stat.h>
#include <dirent.h>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <set>

#include <boost/foreach.hpp>

#include "tools/use-stl-algo-boost-lambda.h"

#include "grid-manager.h"
#include "tools/algorithms.h"
#include "db/abstract-db-connections.h"
#include "tools/read-ini.h"
#include "tools/helper.h"
#include "grid+.h"

#define LOKI_OBJECT_LEVEL_THREADING
#include "loki/Threads.h"

using namespace Grids;
using namespace std;
using namespace Tools;

namespace
{
	struct L : public Loki::ObjectLevelLockable<L> {};
}

VirtualGrid::~VirtualGrid()
{
	for_each(_availableGrids.begin(), _availableGrids.end(),
					 [](GridP* g){ delete g; });//boost::lambda::bind(delete_ptr(), _1));
}

vector<const GridP*> VirtualGrid::availableGrids()
{
	vector<const GridP*> ags(_availableGrids.size());
	std::transform(_availableGrids.begin(), _availableGrids.end(),
								 ags.begin(), boost::lambda::_1);
	return ags;
}

string VirtualGrid::toShortDescription() const
{
	ostringstream s;
	s << "Nutzer-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
	return s.str();
}

bool VirtualGrid::areGridDatasetsAvailable(const list<string>& gdsns) //const
{
	set<string> given;
	std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
								 boost::lambda::_1);

	set<string> all;
  const vector<const GridP*>& ags = availableGrids();
  std::transform(ags.begin(), ags.end(), inserter(all, all.end()),
								 std::bind(&GridP::datasetName, std::placeholders::_1));

	return std::includes(all.begin(), all.end(), given.begin(), given.end());
}

map<string, const GridP*>
    VirtualGrid::gridsForDatasetNames(const list<string>& gdsns) //const
{
	set<string> given;
	std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
								 boost::lambda::_1);

	map<string, const GridP*> m;
  const vector<const GridP*>& ags = availableGrids();
  BOOST_FOREACH(const GridP* ag, ags)
  {
    const string& s = ag->datasetName();
    if(given.find(s) != given.end())
      m.insert(make_pair(s, ag));
	}

	return m;
}

const LatLngPolygonsMatrix& VirtualGrid::latLngCellPolygons()
{
  if(_cellPolygons.isEmpty() && _cols > 0 && _rows > 0)
  {
		//_cellPolygons.resize((_cols-1)*(_rows-1));

		vector<RectCoord> rccs((_cols+1) * (_rows+1));
		//paint the grid cells
    for(unsigned int j = 0; j < _rows+1; j++)
    {
      for(unsigned int i = 0; i < _cols+1; i++)
      {
				double left = _rect.tl.r + (double(i) * _cellSize);
				double top = _rect.tl.h - (double(j) * _cellSize);
        rccs[(j*(_cols+1))+i] = RectCoord(coordinateSystem(), left, top);
				//cout << "(" << left << "," << top << ") ";
			}
		}
		//cout << endl;
		const vector<LatLngCoord>& llcs = RC2latLng(rccs);

		_cellPolygons.resize(_rows, _cols);
    for(unsigned int j = 0; j < _rows; j++)
    {
      for(unsigned int i = 0; i < _cols; i++)
      {
				_cellPolygons[j][i] =
				(LatLngPolygon(llcs.at((j*(_cols+1))+i),
				               llcs.at((j*(_cols+1))+i+1),
				               llcs.at(((j+1)*(_cols+1))+i+1),
				               llcs.at(((j+1)*(_cols+1))+i)));
			}
		}

	}

	return _cellPolygons;
}

//------------------------------------------------------------------------------

map<string, double> VirtualGrid2::dataAt(unsigned int row, unsigned int col) const
{
  map<string, double> res;

  BOOST_FOREACH(Dsn2GridPPtr::value_type p, _dsn2grid)
  {
    res[p.first] = p.second->dataAt(row, col);
  }

  return res;
}

map<string, double> VirtualGrid2::dataAt(const RectCoord& rcc) const
{
  map<string, double> res;

  BOOST_FOREACH(Dsn2GridPPtr::value_type p, _dsn2grid)
  {
    res[p.first] = p.second->dataAt(rcc);
  }

  return res;
}

vector<const GridP*> VirtualGrid2::availableGrids()
{
  vector<const GridP*> ags(_dsn2grid.size());
  std::transform(_dsn2grid.begin(), _dsn2grid.end(),
                 ags.begin(), [](pair<string, GridPPtr> p) { return p.second.get(); });
  return ags;
}

string VirtualGrid2::toShortDescription() const
{
  ostringstream s;
  s << "Nutzer-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
  return s.str();
}

bool VirtualGrid2::areGridDatasetsAvailable(const list<string>& gdsns) //const
{
  set<string> given;
  std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
                 [](string s){return s;});

  set<string> all;
  const vector<const GridP*>& ags = availableGrids();
  std::transform(ags.begin(), ags.end(), inserter(all, all.end()),
                 [](const GridP* g){return g->datasetName();});

  return std::includes(all.begin(), all.end(), given.begin(), given.end());
}

map<string, const GridP*>
VirtualGrid2::gridsForDatasetNames(const list<string>& gdsns) //const
{
  set<string> given;
  std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
                 [](string s){return s;});

  map<string, const GridP*> m;
  const vector<const GridP*>& ags = availableGrids();
  BOOST_FOREACH(const GridP* ag, ags)
  {
    const string& s = ag->datasetName();
    if(given.find(s) != given.end())
      m.insert(make_pair(s, ag));
  }

  return m;
}

const LatLngPolygonsMatrix& VirtualGrid2::latLngCellPolygons()
{
  if(_cellPolygons.isEmpty() && _cols > 0 && _rows > 0)
  {
    //_cellPolygons.resize((_cols-1)*(_rows-1));

    vector<RectCoord> rccs((_cols+1) * (_rows+1));
    //paint the grid cells
    for(unsigned int j = 0; j < _rows+1; j++)
    {
      for(unsigned int i = 0; i < _cols+1; i++)
      {
        double left = _rect.tl.r + (double(i) * _cellSize);
        double top = _rect.tl.h - (double(j) * _cellSize);
        rccs[(j*(_cols+1))+i] = RectCoord(coordinateSystem(), left, top);
        //cout << "(" << left << "," << top << ") ";
      }
    }
    //cout << endl;
    const vector<LatLngCoord>& llcs = RC2latLng(rccs);

    _cellPolygons.resize(_rows, _cols);
    for(unsigned int j = 0; j < _rows; j++)
    {
      for(unsigned int i = 0; i < _cols; i++)
      {
        _cellPolygons[j][i] =
        (LatLngPolygon(llcs.at((j*(_cols+1))+i),
                       llcs.at((j*(_cols+1))+i+1),
                       llcs.at(((j+1)*(_cols+1))+i+1),
                       llcs.at(((j+1)*(_cols+1))+i)));
      }
    }

  }

  return _cellPolygons;
}


//------------------------------------------------------------------------------

NoVirtualGrid::NoVirtualGrid(CoordinateSystem cs,
														 const std::vector<GridProxyPtr>* gps,
														 const Grids::RCRect& rect,
														 double cellSize, unsigned int rows,
														 unsigned int cols,
														 int noDataValue)
	: VirtualGrid(cs, rect, cellSize, rows, cols, noDataValue),
		_gps(gps)
{
}

vector<const GridP*> NoVirtualGrid::availableGrids()
{
  if(_availableGrids.empty())
  {
    BOOST_FOREACH(GridProxyPtr gp, *_gps)
    {
			_availableGrids.push_back(gp->copyOfFullGrid());
		}
	}
	return VirtualGrid::availableGrids();//_availableGrids;
}

string NoVirtualGrid::toShortDescription() const
{
	ostringstream s;
	s << Tools::capitalize(_name) << " ZG:" << cellSize();
	//s << "Standard-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
	return s.str();
}

vector<VirtualGrid::Data> NoVirtualGrid::dataAt(const RectCoord& rcc) const
{
	int ix = int(std::floor((rcc.r - _rect.tl.r) / _cellSize));
	if(ix == int(_cols)) ix--;
	int iy = _rows - int(std::floor((rcc.h - _rect.br.h) / _cellSize)) - 1;
	if(iy < 0) iy = 0;

	return dataAt(iy, ix);
}


//------------------------------------------------------------------------------

NoVirtualGrid2::NoVirtualGrid2(CoordinateSystem cs,
                             const std::vector<GridProxyPtr>* gps,
                             const Grids::RCRect& rect,
                             double cellSize, unsigned int rows,
                             unsigned int cols,
                             int noDataValue)
  : VirtualGrid2(cs, rect, cellSize, rows, cols, map<string, GridPPtr>(), noDataValue),
    _gps(gps)
{
}

vector<const GridP*> NoVirtualGrid2::availableGrids()
{
  if(_dsn2grid.empty())
  {
    BOOST_FOREACH(GridProxyPtr gp, *_gps)
    {
      _dsn2grid[gp->datasetName] = GridPPtr(gp->copyOfFullGrid());
    }
  }
  return VirtualGrid2::availableGrids();//_availableGrids;
}

string NoVirtualGrid2::toShortDescription() const
{
  ostringstream s;
  s << Tools::capitalize(_name) << " ZG:" << cellSize();
  //s << "Standard-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
  return s.str();
}

//------------------------------------------------------------------------------

RealVirtualGrid::RealVirtualGrid(CoordinateSystem cs, const RCRect& rect,
																 double cellSize,
                                 unsigned int rows, unsigned int cols,
																 const VVGP& proxies, int noDataValue)
	: VirtualGrid(cs, rect, cellSize, rows, cols, noDataValue),
		_data(rows, cols),
		_usedGridProxies(proxies.begin(), proxies.end())
{
	//for_each(_data.begin(), _data.end(),
	//         _1 = boost::lambda::bind(new_ptr<vector<Data> >(), _cols));
}

RealVirtualGrid::~RealVirtualGrid()
{
	for_each(_usedGridProxies.begin(), _usedGridProxies.end(),
					 [](vector<GridProxyPtr>* gs){ delete gs; });
					 //boost::lambda::bind(delete_ptr(), _1));
	//for_each(_data.begin(), _data.end(), boost::lambda::bind(delete_ptr(), _1));
}

vector<RealVirtualGrid::Data>
RealVirtualGrid::dataAt(unsigned int row, unsigned int col) const
{
	const vector<Data>& v = _data.valueAt(row, col);
	return v.empty() ? vector<Data>(1, Data()) : v;
}

vector<const GridP*> RealVirtualGrid::availableGrids()
{
  if(_availableGrids.empty())
  {
		map<string, GridP*> dsn2g;

    for(unsigned int i = 0; i < rows(); i++)
    {
      for(unsigned int j = 0; j < cols(); j++)
      {
        BOOST_FOREACH(Data d, dataAt(i, j))
        {
          if(!d.isNoData())
          {
            BOOST_FOREACH(GridProxyPtr gp, *d.gridProxies())
            {
              //target grid
              GridP* tg = valueD(dsn2g, gp->datasetName,
                                 static_cast<GridP*>(NULL));
              if(!tg)
              {
								RectCoord bl = rcRect().toTlTrBrBlVector().at(3);
                tg = new GridP(gp->datasetName, rows(), cols(), cellSize(),
															 bl.r, bl.h, noDataValue(),
															 gp->coordinateSystem);
                dsn2g[gp->datasetName] = tg;
								_availableGrids.push_back(tg);
              }

              tg->setDataAt(i, j, gp->gridPtr()->dataAt(d.row(), d.col()));
						}
					}
				}
			}
		}
	}

	return VirtualGrid::availableGrids();
}

//------------------------------------------------------------------------------

GridManager::GridManager(Env env)
: _env(env)
{
	init("");

	readRegionalizedData();
}

GridManager::~GridManager(){}

void GridManager::init(const Path& userSubPath)
{
	//cout << "entering GridManager::init(" << userSubPath << ")" << endl;

	DIR* dp;
	struct dirent* ep;

	string pathToGrids = _env.asciiGridsPath +
		(userSubPath.empty() ? "" : "/" + userSubPath);
	
	dp = opendir(pathToGrids.c_str());
	if(!dp)
		return;

	//bool gridsInDir = false;
  while((ep = readdir(dp)))
  {
		string dname(ep->d_name);

		//if there are any ascii grids in the top level folder run
		//process there too
		//if(dname.size() > 4 && dname.rfind(".asc") == dname.length() - 4){
		//	gridsInDir = true;
		//} else
		if(dname[0] != '.' && (dname.size() <= 4 ||
													 dname.rfind(".asc") != dname.length() - 4))
		{
			//filter out files starting with . and all files ending in .asc
			//so we basically should get folders
			string folderName(ep->d_name);
			string pathToFolder =
				(userSubPath.empty() ? string("") : userSubPath + "/") + folderName;
			//cout << "diving down into " << pathToFolder << endl;
			init(pathToFolder);
		}
	}

	//cout << "reading grids in " << pathToGrids << endl;
	readGrid2HdfMappingFile(userSubPath);
	checkAndUpdateHdfStore(userSubPath);

	//cout << "leaving GridManager::init(" << userSubPath << ")" << endl;
}

/*
VirtualGrid* GridManager::
createVirtualGrid(const Quadruple<LatLngCoord>& llrect, double cellSize,
									const Path& userSubPath, CoordinateSystem tcs)
{
	const vector<RectCoord>& rcs =
			latLng2RC(asTlTrBrBl<vector<LatLngCoord> >(llrect), tcs);
//	cout << "llrect after conversion to rcs" << endl;
//	for_each(rcs.begin(), rcs.end(), [](RectCoord rc){ cout << rc.toString() << endl; });
	return createVirtualGrid(Quadruple<RectCoord>(rcs), cellSize, userSubPath);
}
*/

//VirtualGrid* GridManager::createVirtualGrid(const Quadruple<RectCoord>& rcpoly,
//                                            double cellSize,
//                                            const Path& userSubPath)
//{
//  Path2GPS::const_iterator ci = _gmdMap.find(userSubPath);
//  if(ci != _gmdMap.end())
//    return createVirtualGrid(ci->second, rcpoly, cellSize);
//  return NULL;
//}

VirtualGrid2* GridManager::createVirtualGrid(const Quadruple<LatLngCoord>& llrect,
                                            const Path& userSubPath)
{
	Path2GPS::const_iterator ci = _gmdMap.find(userSubPath);
	if(ci != _gmdMap.end())
    return createVirtualGrid(ci->second, llrect);//, cellSize);
	return NULL;
}

/*
VirtualGrid*
GridManager::createVirtualGrid(const GMD2GPS& gmd2gridProxies,
															 const Quadruple<RectCoord>& rcpoly,
															 double cellSize)
{
	CoordinateSystem usedCS = rcpoly.tl.coordinateSystem;

	//get bounding rect of rc polygon
	RectCoord tl(usedCS, min(rcpoly.tl.r, rcpoly.bl.r), max(rcpoly.tl.h, rcpoly.tr.h));
	RectCoord br(usedCS, max(rcpoly.tr.r, rcpoly.br.r), min(rcpoly.bl.h, rcpoly.br.h));
	RCRect boundingRect(tl, br);

	vector<GridMetaData> gmds;
	//filter all GridMetaData with correct cellsize and intersecting rect
	BOOST_FOREACH(GMD2GPS::value_type p, gmd2gridProxies)
	{
		//cout << "checking gmd: " << it->first.toString() << endl;
		if(p.first.cellsize == cellSize
			&& p.first.rcRect().intersects(boundingRect))
			gmds.push_back(p.first);
	}

	//cout << "the filtered gmds" << endl;
	//for_each(gmds.begin(), gmds.end(),
	//         cout << boost::lambda::bind(&GridMetaData::toString, _1) << '\n');

	if(gmds.empty())
		return NULL;

	//take the first region of the filtered ones and use this one as
	//base for expanding the bounding rect to full cell-size bounds
	const GridMetaData& firstGmd = gmds.front();
	//cout << "choosen gmd: " << firstGmd.toString() << endl;

	//first find position in choosen grid metadata
	//might in a grid if top left corner of bounding rect is inside
	//the gmd or 0,0 (aka the top left of gmd) if the top left corner of
	//the bounding rect is outside the gmd
	const RCRect& intersectedRect = firstGmd.rcRect().intersected(boundingRect);
	//cout << "intersectedRect: " << intersectedRect.toString() << endl;
	const RectCoord& delta1 = intersectedRect.tl - firstGmd.topLeftCorner();
	int indexR = int(std::floor(delta1.r / cellSize));
	int indexH = int(std::floor(abs(delta1.h / cellSize)));
	//cout << "delta1: " << delta1.toString() << " indexR: " << indexR
	//<< " indexH: " << indexH << endl;
	//rc position into choosen grid-class
	RectCoord firstGmdTl(usedCS,
											 firstGmd.topLeftCorner().r + (indexR * cellSize),
											 firstGmd.topLeftCorner().h - (indexH * cellSize));
	//cout << "grid-class top left: " << firstGmdTl.toString() << endl;

	//now expand the top left corner of the bounding rect to the full
	//outer (hypethetical) grid-cell bound
	//this might be either the same as the previous calculated corner
	//of a grid cell in the gmd (or gmd's 0,0 is case of an exact match)
	//or the bounds have to be extended if the top left corner of bounding rect
	//(of the users selection) was originally outside of the gmd
	RectCoord delta2 = boundingRect.tl - firstGmdTl;
	int nocsToTlr = delta1.r > 0 ? 0 : int(std::ceil(abs(delta2.r / cellSize))); //no of cells
	int nocsToTlh = delta1.h < 0 ? 0 : int(std::ceil(abs(delta2.h / cellSize)));
	//cout << "delta2: " << delta2.toString() << " nocsToTlr: " << nocsToTlr
	//<< " nocsToTlh: " << nocsToTlh << endl;
	//extended tl
	RectCoord etl(usedCS,
								firstGmdTl.r - (nocsToTlr * cellSize),
								firstGmdTl.h + (nocsToTlh * cellSize));
	//cout << "extended top left: " << etl.toString() << endl;

	//adjust br to multiple of cellSize
	RectCoord delta3 = br - etl;
	int nocsR = int(std::ceil(abs(delta3.r / cellSize))); //no of cells
  if(nocsR == 0)
    nocsR++; //the selection is choosing at least one cell
	int nocsH = int(std::ceil(abs(delta3.h / cellSize)));
  if(nocsH == 0)
    nocsH++;
	//cout << "delta3: " << delta3.toString() << " nocsR: " << nocsR
	//<< " nocsH: " << nocsH << endl;
	//the extended final bounding rect
	RectCoord ebr(usedCS,
								etl.r + (nocsR * cellSize),
								etl.h - (nocsH * cellSize));
	RCRect extendedBoundingRect(etl, ebr);
	//cout << "extended bounding rect: " << extendedBoundingRect.toString() << endl;

	//for efficiency reasons every data element just references a vector
	//of gridproxies with available grids from a given region (same gridmetadata)
	//but these have to be unique, thus a number of vectors is created below
	//which will be stored in the according virtual grid which can
	//be referenced directly in the data elements
	map<GridMetaData, GridProxies*> gmd2gps;
	//get the unique set of all dataset names to be used below
	set<string> uniqueDatasetNames;
  BOOST_FOREACH(const GridMetaData& gmd, gmds)
  {
		GridProxies* gps = new GridProxies;
		GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
		if(ci == gmd2gridProxies.end())
			continue;
		const GridProxies& agps = ci->second;
    BOOST_FOREACH(GridProxyPtr agp, agps)
    {
      string s = agp->datasetName;
      if(uniqueDatasetNames.find(s) == uniqueDatasetNames.end())
      {
				uniqueDatasetNames.insert(s);
        gps->push_back(agp);
			}
		}

		gmd2gps.insert(make_pair(gmd, gps));

		//cout << "copied these proxies: " << endl;
		//for_each(gps->begin(), gps->end(),
		//         cout << boost::lambda::bind(&GridProxy::toString, _1) << "\n");
	}
	//for_each(uniqueDatasetNames.begin(), uniqueDatasetNames.end(),
	//         cout << _1 << '\n');
	vector<GridProxies*> gpss;
	transform(gmd2gps.begin(), gmd2gps.end(), back_inserter(gpss),
						[](pair<const GridMetaData, GridProxies*> v){ return v.second; }
						//boost::lambda::bind<GridProxies*>(&pair<const GridMetaData, GridProxies*>::second,
															 //_1)
																	 );

	RealVirtualGrid* vg = new RealVirtualGrid(usedCS, extendedBoundingRect,
																						cellSize,
	                                          nocsH, nocsR, gpss);

	//now we have to iterate through the virtual grid and fill its cells
	//either with no data values or with references to the potential grids
	//and the correct indices into it
  for(unsigned int i = 0; i < vg->rows(); i++)
  {
    for(unsigned int k = 0; k < vg->cols(); k++)
    {
			vector<GridMetaData>::iterator it = gmds.begin();// - 1;
			const RectCoord& c = vg->rcCoordAt(i, k);

			auto f = [&](const GridMetaData& gmd){ return gmd.rcRect().contains(c, true); };

			while((it = find_if(it, gmds.end(), f)) != gmds.end())
			{
				//else it vg is preinitialized with no data values
				const GridMetaData& gmd = *it;

				//this is just temporary as it won't work if a virtual grid
				//consists of data from more than one region, but better than
				//nothing until the whole virtual grid thing will be redone
				if(gmd.regionName == "uecker" ||
					 gmd.regionName == "sachsen" ||
					 gmd.regionName == "brazil-sinop" ||
					 gmd.regionName == "brazil-campo-verde")
					vg->setCustomId(gmd.regionName);

				pair<Row, Col> rc = rowColInGrid(gmd, c);
				//cout << "using gmd: " << gmd.toString() << endl;
				//cout << "( " << i << "," << k << ") indexR (col): " << rc.col
				//<< " indexH (row): " << rc.row << endl;

				//cout << "adding data at: " << i << "/" << k << endl;
        vg->addDataAt(i, k, VirtualGrid::Data(gmd2gps[gmd],
                                              rc.first, rc.second));

				it++;
			}
		}
	}

	return vg;
}
*/

VirtualGrid2*
GridManager::createVirtualGrid(const GMD2GPS& gmd2gridProxies,
                               const Quadruple<Tools::LatLngCoord>& llrect)
{
  //map<CoordinateSystem, RCRect> cs2boundingRect;
  typedef map<CoordinateSystem, vector<GridMetaData>> CS2GMDS;
  CS2GMDS cs2gmds;
  int minCellSize = 100000000;
  typedef map<CoordinateSystem, int> CS2NO;
  CS2NO cs2noOfGPS;
  //filter all GridMetaData with correct cellsize and intersecting rect
  BOOST_FOREACH(GMD2GPS::value_type p, gmd2gridProxies)
  {
    auto cs = p.first.coordinateSystem;
//    auto ci = cs2boundingRect.find(cs);
//    RCRect boundingRect;
//    if(ci == cs2boundingRect.end())
//    {
      //first convert latlng to rect coordinates of gmd
      const vector<RectCoord>& rcs = latLng2RC(asTlTrBrBl<vector<LatLngCoord> >(llrect), cs);
      auto rcpoly = Quadruple<RectCoord>(rcs);

      //create and store a bounding rect in this rect coordinate system
      RectCoord tl(cs, min(rcpoly.tl.r, rcpoly.bl.r), max(rcpoly.tl.h, rcpoly.tr.h));
      RectCoord br(cs, max(rcpoly.tr.r, rcpoly.br.r), min(rcpoly.bl.h, rcpoly.br.h));
      RCRect boundingRect = RCRect(tl, br);
//      cs2boundingRect[cs] = boundingRect;
//    }
//    else
//      boundingRect = ci->second;

    //cout << "checking gmd: " << it->first.toString() << endl;
    if(p.first.rcRect().intersects(boundingRect))
    {
      cs2noOfGPS[cs] += p.second.size();
      minCellSize = min(p.first.cellsize, minCellSize);
      cs2gmds[cs].push_back(p.first);
    }
  }

  //find most often used coordinate system in the selected region
  //we assume that this coordinate system can safely be used in this region,
  //else there wouldn't be that many grids there
  CoordinateSystem usedCS = CoordinateSystem();
  int noOfGPS = 0;
  BOOST_FOREACH(CS2NO::value_type p, cs2noOfGPS)
  {
    if(p.second > noOfGPS)
    {
      noOfGPS = p.second;
      usedCS = p.first;
    }
  }

  //if we didn't find any coordinate system we leave
  if(usedCS == CoordinateSystem())
    return NULL;

  //determine the dimensions of the bounding rect if the top left corner is exactly
  //the user choosen one and the bottom right corner is the one determined to include the user choosen one
  //given the determined cellsize
  const vector<RectCoord>& rcs = latLng2RC(asTlTrBrBl<vector<LatLngCoord> >(llrect), usedCS);
  auto rcpoly = Quadruple<RectCoord>(rcs);
  int noOfRows = int(std::ceil((rcpoly.tl.h - rcpoly.br.h)/minCellSize));
  int noOfCols = int(std::ceil((rcpoly.br.r - rcpoly.tl.r)/minCellSize));
  RectCoord tl(usedCS, min(rcpoly.tl.r, rcpoly.bl.r), max(rcpoly.tl.h, rcpoly.tr.h));
  RectCoord br(usedCS, tl.r + noOfCols*minCellSize, tl.h - noOfRows*minCellSize);

  map<string, GridPPtr> dsn2grid;
  set<string> uniqueDatasetNames;
  BOOST_FOREACH(CS2GMDS::value_type p, cs2gmds)
  {
    BOOST_FOREACH(const GridMetaData& gmd, p.second)
    {
      GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
      if(ci == gmd2gridProxies.end())
        continue;
      const GridProxies& agps = ci->second;
      BOOST_FOREACH(GridProxyPtr agp, agps)
      {
//        cout << "dsn: |" << agp->datasetName << "|" << endl;
        uniqueDatasetNames.insert(agp->datasetName);
      }
    }
  }

  //create all the grids needed
  BOOST_FOREACH(string dsn, uniqueDatasetNames)
  {
    dsn2grid[dsn] = GridPPtr(new GridP(dsn, noOfRows, noOfCols, minCellSize, tl.r, br.h, -9999, usedCS));
  }

  if(dsn2grid.empty())
    return NULL;

  string regionName;

  GridPPtr someGrid = dsn2grid.begin()->second;
  for(int r = 0, rs = someGrid->rows(); r < rs; r++)
  {
    for(int c = 0, cs = someGrid->cols(); c < cs; c++)
    {
      RectCoord cellCenter = someGrid->rcCoordAtCenter(r, c);
      double sgCellSize = someGrid->cellSize();

      LatLngCoord llcc = RC2latLng(cellCenter);

      BOOST_FOREACH(CS2GMDS::value_type p, cs2gmds)
      {
        CoordinateSystem cs = p.first;
        RectCoord rccc = latLng2RC(llcc, cs);

        BOOST_FOREACH(const GridMetaData& gmd, p.second)
        {
//          cout << "gmd: " << gmd.toCanonicalString() << " regionName: " << regionName << endl;
          regionName = gmd.regionName;
//          if(gmd.regionName == "uecker" ||
//             gmd.regionName == "sachsen" ||
//             gmd.regionName == "uelzen" ||
//             gmd.regionName == "brazil-sinop" ||
//             gmd.regionName == "brazil-campo-verde")
//            regionName = gmd.regionName;

          GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
          if(ci == gmd2gridProxies.end())
            continue;
          const GridProxies& agps = ci->second;
          if(agps.empty())
            continue;

          GridP* g = agps.front()->gridPtr();
          double gCellSize = g->cellSize();

          auto pos = g->rc2rowCol(rccc);

          //if any value lies farther outside than a cellwidth of the target grid, it can't intersect the closest
          //source cell thus we associate a no data value with the target grid value
          bool setNoDataValue =
              (pos.row < 0 && abs(pos.row*gCellSize) > sgCellSize) ||
              (pos.col < 0 && abs(pos.col*gCellSize) > sgCellSize) ||
              (pos.row > 0 && pos.isRowOutside && abs((pos.row - g->rows())*gCellSize > sgCellSize)) ||
              (pos.col > 0 && pos.isColOutside && abs((pos.col - g->cols())*gCellSize > sgCellSize));

          //if the position is in the cell outside, just associate the target grid value with the border value of
          //the source grid
          if(!setNoDataValue)
          {
            if(pos.row < 0)
              pos.row = 0;
            else if(pos.row > 0 && pos.isRowOutside)
              pos.row = g->rows() - 1;
            if(pos.col < 0)
              pos.col = 0;
            else if(pos.col > 0 && pos.isColOutside)
              pos.col = g->cols() - 1;
          }

          BOOST_FOREACH(GridProxyPtr agp, agps)
          {
            GridPPtr tg = dsn2grid[agp->datasetName];
            if(setNoDataValue)
              tg->setNoDataValueAt(r,c);
            else
              tg->setDataAt(r, c, agp->gridPtr()->dataAt(pos.row, pos.col));
          }
        }
      }
    }
  }

  auto vg = new VirtualGrid2(usedCS, RCRect(tl, br), minCellSize, noOfRows, noOfCols, dsn2grid);
  vg->setCustomId(regionName);

  return vg;
}

VirtualGrid2*
GridManager::virtualGridForGridMetaData(const GridMetaData& gmd,
                                        const Path& userSubPath)
{
	Path2GPS::const_iterator ci = _gmdMap.find(userSubPath);
  if(ci != _gmdMap.end())
  {
		GMD2GPS::const_iterator ci2 = ci->second.find(gmd);
    if(ci2 != ci->second.end())
    {
      NoVirtualGrid2* nvg =
          new NoVirtualGrid2(gmd.coordinateSystem,
														&(ci2->second),
														gmd.rcRect(), gmd.cellsize,
                            gmd.nrows, gmd.ncols);
			//cout << "gmd.regionName: " << gmd.regionName << endl;
			nvg->setName(gmd.regionName);
			nvg->setCustomId(gmd.regionName);
			return nvg;
		}
	}<commit_msg>test to speedup creation of virtual grids<commit_after>/**
Authors:
<NAME> <<EMAIL>>

Maintainers:
Currently maintained by the authors.

This file is part of the util library used by models created at the Institute of
Landscape Systems Analysis at the ZALF.
Copyright (C) 2007-2013, Leibniz Centre for Agricultural Landscape Research (ZALF)

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <fstream>
#include <sstream>
#include <cstddef>
#include <cstdio>
#include <sys/types.h>
#include <sys/stat.h>
#include <dirent.h>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <set>

#include <boost/foreach.hpp>

#include "tools/use-stl-algo-boost-lambda.h"

#include "grid-manager.h"
#include "tools/algorithms.h"
#include "db/abstract-db-connections.h"
#include "tools/read-ini.h"
#include "tools/helper.h"
#include "grid+.h"

#define LOKI_OBJECT_LEVEL_THREADING
#include "loki/Threads.h"

using namespace Grids;
using namespace std;
using namespace Tools;

namespace
{
	struct L : public Loki::ObjectLevelLockable<L> {};
}

VirtualGrid::~VirtualGrid()
{
	for_each(_availableGrids.begin(), _availableGrids.end(),
					 [](GridP* g){ delete g; });//boost::lambda::bind(delete_ptr(), _1));
}

vector<const GridP*> VirtualGrid::availableGrids()
{
	vector<const GridP*> ags(_availableGrids.size());
	std::transform(_availableGrids.begin(), _availableGrids.end(),
								 ags.begin(), boost::lambda::_1);
	return ags;
}

string VirtualGrid::toShortDescription() const
{
	ostringstream s;
	s << "Nutzer-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
	return s.str();
}

bool VirtualGrid::areGridDatasetsAvailable(const list<string>& gdsns) //const
{
	set<string> given;
	std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
								 boost::lambda::_1);

	set<string> all;
  const vector<const GridP*>& ags = availableGrids();
  std::transform(ags.begin(), ags.end(), inserter(all, all.end()),
								 std::bind(&GridP::datasetName, std::placeholders::_1));

	return std::includes(all.begin(), all.end(), given.begin(), given.end());
}

map<string, const GridP*>
    VirtualGrid::gridsForDatasetNames(const list<string>& gdsns) //const
{
	set<string> given;
	std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
								 boost::lambda::_1);

	map<string, const GridP*> m;
  const vector<const GridP*>& ags = availableGrids();
  BOOST_FOREACH(const GridP* ag, ags)
  {
    const string& s = ag->datasetName();
    if(given.find(s) != given.end())
      m.insert(make_pair(s, ag));
	}

	return m;
}

const LatLngPolygonsMatrix& VirtualGrid::latLngCellPolygons()
{
  if(_cellPolygons.isEmpty() && _cols > 0 && _rows > 0)
  {
		//_cellPolygons.resize((_cols-1)*(_rows-1));

		vector<RectCoord> rccs((_cols+1) * (_rows+1));
		//paint the grid cells
    for(unsigned int j = 0; j < _rows+1; j++)
    {
      for(unsigned int i = 0; i < _cols+1; i++)
      {
				double left = _rect.tl.r + (double(i) * _cellSize);
				double top = _rect.tl.h - (double(j) * _cellSize);
        rccs[(j*(_cols+1))+i] = RectCoord(coordinateSystem(), left, top);
				//cout << "(" << left << "," << top << ") ";
			}
		}
		//cout << endl;
		const vector<LatLngCoord>& llcs = RC2latLng(rccs);

		_cellPolygons.resize(_rows, _cols);
    for(unsigned int j = 0; j < _rows; j++)
    {
      for(unsigned int i = 0; i < _cols; i++)
      {
				_cellPolygons[j][i] =
				(LatLngPolygon(llcs.at((j*(_cols+1))+i),
				               llcs.at((j*(_cols+1))+i+1),
				               llcs.at(((j+1)*(_cols+1))+i+1),
				               llcs.at(((j+1)*(_cols+1))+i)));
			}
		}

	}

	return _cellPolygons;
}

//------------------------------------------------------------------------------

map<string, double> VirtualGrid2::dataAt(unsigned int row, unsigned int col) const
{
  map<string, double> res;

  BOOST_FOREACH(Dsn2GridPPtr::value_type p, _dsn2grid)
  {
    res[p.first] = p.second->dataAt(row, col);
  }

  return res;
}

map<string, double> VirtualGrid2::dataAt(const RectCoord& rcc) const
{
  map<string, double> res;

  BOOST_FOREACH(Dsn2GridPPtr::value_type p, _dsn2grid)
  {
    res[p.first] = p.second->dataAt(rcc);
  }

  return res;
}

vector<const GridP*> VirtualGrid2::availableGrids()
{
  vector<const GridP*> ags(_dsn2grid.size());
  std::transform(_dsn2grid.begin(), _dsn2grid.end(),
                 ags.begin(), [](pair<string, GridPPtr> p) { return p.second.get(); });
  return ags;
}

string VirtualGrid2::toShortDescription() const
{
  ostringstream s;
  s << "Nutzer-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
  return s.str();
}

bool VirtualGrid2::areGridDatasetsAvailable(const list<string>& gdsns) //const
{
  set<string> given;
  std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
                 [](string s){return s;});

  set<string> all;
  const vector<const GridP*>& ags = availableGrids();
  std::transform(ags.begin(), ags.end(), inserter(all, all.end()),
                 [](const GridP* g){return g->datasetName();});

  return std::includes(all.begin(), all.end(), given.begin(), given.end());
}

map<string, const GridP*>
VirtualGrid2::gridsForDatasetNames(const list<string>& gdsns) //const
{
  set<string> given;
  std::transform(gdsns.begin(), gdsns.end(), inserter(given, given.end()),
                 [](string s){return s;});

  map<string, const GridP*> m;
  const vector<const GridP*>& ags = availableGrids();
  BOOST_FOREACH(const GridP* ag, ags)
  {
    const string& s = ag->datasetName();
    if(given.find(s) != given.end())
      m.insert(make_pair(s, ag));
  }

  return m;
}

const LatLngPolygonsMatrix& VirtualGrid2::latLngCellPolygons()
{
  if(_cellPolygons.isEmpty() && _cols > 0 && _rows > 0)
  {
    //_cellPolygons.resize((_cols-1)*(_rows-1));

    vector<RectCoord> rccs((_cols+1) * (_rows+1));
    //paint the grid cells
    for(unsigned int j = 0; j < _rows+1; j++)
    {
      for(unsigned int i = 0; i < _cols+1; i++)
      {
        double left = _rect.tl.r + (double(i) * _cellSize);
        double top = _rect.tl.h - (double(j) * _cellSize);
        rccs[(j*(_cols+1))+i] = RectCoord(coordinateSystem(), left, top);
        //cout << "(" << left << "," << top << ") ";
      }
    }
    //cout << endl;
    const vector<LatLngCoord>& llcs = RC2latLng(rccs);

    _cellPolygons.resize(_rows, _cols);
    for(unsigned int j = 0; j < _rows; j++)
    {
      for(unsigned int i = 0; i < _cols; i++)
      {
        _cellPolygons[j][i] =
        (LatLngPolygon(llcs.at((j*(_cols+1))+i),
                       llcs.at((j*(_cols+1))+i+1),
                       llcs.at(((j+1)*(_cols+1))+i+1),
                       llcs.at(((j+1)*(_cols+1))+i)));
      }
    }

  }

  return _cellPolygons;
}


//------------------------------------------------------------------------------

NoVirtualGrid::NoVirtualGrid(CoordinateSystem cs,
														 const std::vector<GridProxyPtr>* gps,
														 const Grids::RCRect& rect,
														 double cellSize, unsigned int rows,
														 unsigned int cols,
														 int noDataValue)
	: VirtualGrid(cs, rect, cellSize, rows, cols, noDataValue),
		_gps(gps)
{
}

vector<const GridP*> NoVirtualGrid::availableGrids()
{
  if(_availableGrids.empty())
  {
    BOOST_FOREACH(GridProxyPtr gp, *_gps)
    {
			_availableGrids.push_back(gp->copyOfFullGrid());
		}
	}
	return VirtualGrid::availableGrids();//_availableGrids;
}

string NoVirtualGrid::toShortDescription() const
{
	ostringstream s;
	s << Tools::capitalize(_name) << " ZG:" << cellSize();
	//s << "Standard-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
	return s.str();
}

vector<VirtualGrid::Data> NoVirtualGrid::dataAt(const RectCoord& rcc) const
{
	int ix = int(std::floor((rcc.r - _rect.tl.r) / _cellSize));
	if(ix == int(_cols)) ix--;
	int iy = _rows - int(std::floor((rcc.h - _rect.br.h) / _cellSize)) - 1;
	if(iy < 0) iy = 0;

	return dataAt(iy, ix);
}


//------------------------------------------------------------------------------

NoVirtualGrid2::NoVirtualGrid2(CoordinateSystem cs,
                             const std::vector<GridProxyPtr>* gps,
                             const Grids::RCRect& rect,
                             double cellSize, unsigned int rows,
                             unsigned int cols,
                             int noDataValue)
  : VirtualGrid2(cs, rect, cellSize, rows, cols, map<string, GridPPtr>(), noDataValue),
    _gps(gps)
{
}

vector<const GridP*> NoVirtualGrid2::availableGrids()
{
  if(_dsn2grid.empty())
  {
    BOOST_FOREACH(GridProxyPtr gp, *_gps)
    {
      _dsn2grid[gp->datasetName] = GridPPtr(gp->copyOfFullGrid());
    }
  }
  return VirtualGrid2::availableGrids();//_availableGrids;
}

string NoVirtualGrid2::toShortDescription() const
{
  ostringstream s;
  s << Tools::capitalize(_name) << " ZG:" << cellSize();
  //s << "Standard-Region [" << cols() << "x" << rows() << "], ZG:" << cellSize();
  return s.str();
}

//------------------------------------------------------------------------------

RealVirtualGrid::RealVirtualGrid(CoordinateSystem cs, const RCRect& rect,
																 double cellSize,
                                 unsigned int rows, unsigned int cols,
																 const VVGP& proxies, int noDataValue)
	: VirtualGrid(cs, rect, cellSize, rows, cols, noDataValue),
		_data(rows, cols),
		_usedGridProxies(proxies.begin(), proxies.end())
{
	//for_each(_data.begin(), _data.end(),
	//         _1 = boost::lambda::bind(new_ptr<vector<Data> >(), _cols));
}

RealVirtualGrid::~RealVirtualGrid()
{
	for_each(_usedGridProxies.begin(), _usedGridProxies.end(),
					 [](vector<GridProxyPtr>* gs){ delete gs; });
					 //boost::lambda::bind(delete_ptr(), _1));
	//for_each(_data.begin(), _data.end(), boost::lambda::bind(delete_ptr(), _1));
}

vector<RealVirtualGrid::Data>
RealVirtualGrid::dataAt(unsigned int row, unsigned int col) const
{
	const vector<Data>& v = _data.valueAt(row, col);
	return v.empty() ? vector<Data>(1, Data()) : v;
}

vector<const GridP*> RealVirtualGrid::availableGrids()
{
  if(_availableGrids.empty())
  {
		map<string, GridP*> dsn2g;

    for(unsigned int i = 0; i < rows(); i++)
    {
      for(unsigned int j = 0; j < cols(); j++)
      {
        BOOST_FOREACH(Data d, dataAt(i, j))
        {
          if(!d.isNoData())
          {
            BOOST_FOREACH(GridProxyPtr gp, *d.gridProxies())
            {
              //target grid
              GridP* tg = valueD(dsn2g, gp->datasetName,
                                 static_cast<GridP*>(NULL));
              if(!tg)
              {
								RectCoord bl = rcRect().toTlTrBrBlVector().at(3);
                tg = new GridP(gp->datasetName, rows(), cols(), cellSize(),
															 bl.r, bl.h, noDataValue(),
															 gp->coordinateSystem);
                dsn2g[gp->datasetName] = tg;
								_availableGrids.push_back(tg);
              }

              tg->setDataAt(i, j, gp->gridPtr()->dataAt(d.row(), d.col()));
						}
					}
				}
			}
		}
	}

	return VirtualGrid::availableGrids();
}

//------------------------------------------------------------------------------

GridManager::GridManager(Env env)
: _env(env)
{
	init("");

	readRegionalizedData();
}

GridManager::~GridManager(){}

void GridManager::init(const Path& userSubPath)
{
	//cout << "entering GridManager::init(" << userSubPath << ")" << endl;

	DIR* dp;
	struct dirent* ep;

	string pathToGrids = _env.asciiGridsPath +
		(userSubPath.empty() ? "" : "/" + userSubPath);
	
	dp = opendir(pathToGrids.c_str());
	if(!dp)
		return;

	//bool gridsInDir = false;
  while((ep = readdir(dp)))
  {
		string dname(ep->d_name);

		//if there are any ascii grids in the top level folder run
		//process there too
		//if(dname.size() > 4 && dname.rfind(".asc") == dname.length() - 4){
		//	gridsInDir = true;
		//} else
		if(dname[0] != '.' && (dname.size() <= 4 ||
													 dname.rfind(".asc") != dname.length() - 4))
		{
			//filter out files starting with . and all files ending in .asc
			//so we basically should get folders
			string folderName(ep->d_name);
			string pathToFolder =
				(userSubPath.empty() ? string("") : userSubPath + "/") + folderName;
			//cout << "diving down into " << pathToFolder << endl;
			init(pathToFolder);
		}
	}

	//cout << "reading grids in " << pathToGrids << endl;
	readGrid2HdfMappingFile(userSubPath);
	checkAndUpdateHdfStore(userSubPath);

	//cout << "leaving GridManager::init(" << userSubPath << ")" << endl;
}

/*
VirtualGrid* GridManager::
createVirtualGrid(const Quadruple<LatLngCoord>& llrect, double cellSize,
									const Path& userSubPath, CoordinateSystem tcs)
{
	const vector<RectCoord>& rcs =
			latLng2RC(asTlTrBrBl<vector<LatLngCoord> >(llrect), tcs);
//	cout << "llrect after conversion to rcs" << endl;
//	for_each(rcs.begin(), rcs.end(), [](RectCoord rc){ cout << rc.toString() << endl; });
	return createVirtualGrid(Quadruple<RectCoord>(rcs), cellSize, userSubPath);
}
*/

//VirtualGrid* GridManager::createVirtualGrid(const Quadruple<RectCoord>& rcpoly,
//                                            double cellSize,
//                                            const Path& userSubPath)
//{
//  Path2GPS::const_iterator ci = _gmdMap.find(userSubPath);
//  if(ci != _gmdMap.end())
//    return createVirtualGrid(ci->second, rcpoly, cellSize);
//  return NULL;
//}

VirtualGrid2* GridManager::createVirtualGrid(const Quadruple<LatLngCoord>& llrect,
                                            const Path& userSubPath)
{
	Path2GPS::const_iterator ci = _gmdMap.find(userSubPath);
	if(ci != _gmdMap.end())
    return createVirtualGrid(ci->second, llrect);//, cellSize);
	return NULL;
}

/*
VirtualGrid*
GridManager::createVirtualGrid(const GMD2GPS& gmd2gridProxies,
															 const Quadruple<RectCoord>& rcpoly,
															 double cellSize)
{
	CoordinateSystem usedCS = rcpoly.tl.coordinateSystem;

	//get bounding rect of rc polygon
	RectCoord tl(usedCS, min(rcpoly.tl.r, rcpoly.bl.r), max(rcpoly.tl.h, rcpoly.tr.h));
	RectCoord br(usedCS, max(rcpoly.tr.r, rcpoly.br.r), min(rcpoly.bl.h, rcpoly.br.h));
	RCRect boundingRect(tl, br);

	vector<GridMetaData> gmds;
	//filter all GridMetaData with correct cellsize and intersecting rect
	BOOST_FOREACH(GMD2GPS::value_type p, gmd2gridProxies)
	{
		//cout << "checking gmd: " << it->first.toString() << endl;
		if(p.first.cellsize == cellSize
			&& p.first.rcRect().intersects(boundingRect))
			gmds.push_back(p.first);
	}

	//cout << "the filtered gmds" << endl;
	//for_each(gmds.begin(), gmds.end(),
	//         cout << boost::lambda::bind(&GridMetaData::toString, _1) << '\n');

	if(gmds.empty())
		return NULL;

	//take the first region of the filtered ones and use this one as
	//base for expanding the bounding rect to full cell-size bounds
	const GridMetaData& firstGmd = gmds.front();
	//cout << "choosen gmd: " << firstGmd.toString() << endl;

	//first find position in choosen grid metadata
	//might in a grid if top left corner of bounding rect is inside
	//the gmd or 0,0 (aka the top left of gmd) if the top left corner of
	//the bounding rect is outside the gmd
	const RCRect& intersectedRect = firstGmd.rcRect().intersected(boundingRect);
	//cout << "intersectedRect: " << intersectedRect.toString() << endl;
	const RectCoord& delta1 = intersectedRect.tl - firstGmd.topLeftCorner();
	int indexR = int(std::floor(delta1.r / cellSize));
	int indexH = int(std::floor(abs(delta1.h / cellSize)));
	//cout << "delta1: " << delta1.toString() << " indexR: " << indexR
	//<< " indexH: " << indexH << endl;
	//rc position into choosen grid-class
	RectCoord firstGmdTl(usedCS,
											 firstGmd.topLeftCorner().r + (indexR * cellSize),
											 firstGmd.topLeftCorner().h - (indexH * cellSize));
	//cout << "grid-class top left: " << firstGmdTl.toString() << endl;

	//now expand the top left corner of the bounding rect to the full
	//outer (hypethetical) grid-cell bound
	//this might be either the same as the previous calculated corner
	//of a grid cell in the gmd (or gmd's 0,0 is case of an exact match)
	//or the bounds have to be extended if the top left corner of bounding rect
	//(of the users selection) was originally outside of the gmd
	RectCoord delta2 = boundingRect.tl - firstGmdTl;
	int nocsToTlr = delta1.r > 0 ? 0 : int(std::ceil(abs(delta2.r / cellSize))); //no of cells
	int nocsToTlh = delta1.h < 0 ? 0 : int(std::ceil(abs(delta2.h / cellSize)));
	//cout << "delta2: " << delta2.toString() << " nocsToTlr: " << nocsToTlr
	//<< " nocsToTlh: " << nocsToTlh << endl;
	//extended tl
	RectCoord etl(usedCS,
								firstGmdTl.r - (nocsToTlr * cellSize),
								firstGmdTl.h + (nocsToTlh * cellSize));
	//cout << "extended top left: " << etl.toString() << endl;

	//adjust br to multiple of cellSize
	RectCoord delta3 = br - etl;
	int nocsR = int(std::ceil(abs(delta3.r / cellSize))); //no of cells
  if(nocsR == 0)
    nocsR++; //the selection is choosing at least one cell
	int nocsH = int(std::ceil(abs(delta3.h / cellSize)));
  if(nocsH == 0)
    nocsH++;
	//cout << "delta3: " << delta3.toString() << " nocsR: " << nocsR
	//<< " nocsH: " << nocsH << endl;
	//the extended final bounding rect
	RectCoord ebr(usedCS,
								etl.r + (nocsR * cellSize),
								etl.h - (nocsH * cellSize));
	RCRect extendedBoundingRect(etl, ebr);
	//cout << "extended bounding rect: " << extendedBoundingRect.toString() << endl;

	//for efficiency reasons every data element just references a vector
	//of gridproxies with available grids from a given region (same gridmetadata)
	//but these have to be unique, thus a number of vectors is created below
	//which will be stored in the according virtual grid which can
	//be referenced directly in the data elements
	map<GridMetaData, GridProxies*> gmd2gps;
	//get the unique set of all dataset names to be used below
	set<string> uniqueDatasetNames;
  BOOST_FOREACH(const GridMetaData& gmd, gmds)
  {
		GridProxies* gps = new GridProxies;
		GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
		if(ci == gmd2gridProxies.end())
			continue;
		const GridProxies& agps = ci->second;
    BOOST_FOREACH(GridProxyPtr agp, agps)
    {
      string s = agp->datasetName;
      if(uniqueDatasetNames.find(s) == uniqueDatasetNames.end())
      {
				uniqueDatasetNames.insert(s);
        gps->push_back(agp);
			}
		}

		gmd2gps.insert(make_pair(gmd, gps));

		//cout << "copied these proxies: " << endl;
		//for_each(gps->begin(), gps->end(),
		//         cout << boost::lambda::bind(&GridProxy::toString, _1) << "\n");
	}
	//for_each(uniqueDatasetNames.begin(), uniqueDatasetNames.end(),
	//         cout << _1 << '\n');
	vector<GridProxies*> gpss;
	transform(gmd2gps.begin(), gmd2gps.end(), back_inserter(gpss),
						[](pair<const GridMetaData, GridProxies*> v){ return v.second; }
						//boost::lambda::bind<GridProxies*>(&pair<const GridMetaData, GridProxies*>::second,
															 //_1)
																	 );

	RealVirtualGrid* vg = new RealVirtualGrid(usedCS, extendedBoundingRect,
																						cellSize,
	                                          nocsH, nocsR, gpss);

	//now we have to iterate through the virtual grid and fill its cells
	//either with no data values or with references to the potential grids
	//and the correct indices into it
  for(unsigned int i = 0; i < vg->rows(); i++)
  {
    for(unsigned int k = 0; k < vg->cols(); k++)
    {
			vector<GridMetaData>::iterator it = gmds.begin();// - 1;
			const RectCoord& c = vg->rcCoordAt(i, k);

			auto f = [&](const GridMetaData& gmd){ return gmd.rcRect().contains(c, true); };

			while((it = find_if(it, gmds.end(), f)) != gmds.end())
			{
				//else it vg is preinitialized with no data values
				const GridMetaData& gmd = *it;

				//this is just temporary as it won't work if a virtual grid
				//consists of data from more than one region, but better than
				//nothing until the whole virtual grid thing will be redone
				if(gmd.regionName == "uecker" ||
					 gmd.regionName == "sachsen" ||
					 gmd.regionName == "brazil-sinop" ||
					 gmd.regionName == "brazil-campo-verde")
					vg->setCustomId(gmd.regionName);

				pair<Row, Col> rc = rowColInGrid(gmd, c);
				//cout << "using gmd: " << gmd.toString() << endl;
				//cout << "( " << i << "," << k << ") indexR (col): " << rc.col
				//<< " indexH (row): " << rc.row << endl;

				//cout << "adding data at: " << i << "/" << k << endl;
        vg->addDataAt(i, k, VirtualGrid::Data(gmd2gps[gmd],
                                              rc.first, rc.second));

				it++;
			}
		}
	}

	return vg;
}
*/

VirtualGrid2*
GridManager::createVirtualGrid(const GMD2GPS& gmd2gridProxies,
                               const Quadruple<Tools::LatLngCoord>& llrect)
{
  //map<CoordinateSystem, RCRect> cs2boundingRect;
  typedef map<CoordinateSystem, vector<GridMetaData>> CS2GMDS;
  CS2GMDS cs2gmds;
  int minCellSize = 100000000;
  typedef map<CoordinateSystem, int> CS2NO;
  CS2NO cs2noOfGPS;
  //filter all GridMetaData with correct cellsize and intersecting rect
  BOOST_FOREACH(GMD2GPS::value_type p, gmd2gridProxies)
  {
    auto cs = p.first.coordinateSystem;
//    auto ci = cs2boundingRect.find(cs);
//    RCRect boundingRect;
//    if(ci == cs2boundingRect.end())
//    {
      //first convert latlng to rect coordinates of gmd
      const vector<RectCoord>& rcs = latLng2RC(asTlTrBrBl<vector<LatLngCoord> >(llrect), cs);
      auto rcpoly = Quadruple<RectCoord>(rcs);

      //create and store a bounding rect in this rect coordinate system
      RectCoord tl(cs, min(rcpoly.tl.r, rcpoly.bl.r), max(rcpoly.tl.h, rcpoly.tr.h));
      RectCoord br(cs, max(rcpoly.tr.r, rcpoly.br.r), min(rcpoly.bl.h, rcpoly.br.h));
      RCRect boundingRect = RCRect(tl, br);
//      cs2boundingRect[cs] = boundingRect;
//    }
//    else
//      boundingRect = ci->second;

    //cout << "checking gmd: " << it->first.toString() << endl;
    if(p.first.rcRect().intersects(boundingRect))
    {
      cs2noOfGPS[cs] += p.second.size();
      minCellSize = min(p.first.cellsize, minCellSize);
      cs2gmds[cs].push_back(p.first);
    }
  }

  //find most often used coordinate system in the selected region
  //we assume that this coordinate system can safely be used in this region,
  //else there wouldn't be that many grids there
  CoordinateSystem usedCS = CoordinateSystem();
  int noOfGPS = 0;
  BOOST_FOREACH(CS2NO::value_type p, cs2noOfGPS)
  {
    if(p.second > noOfGPS)
    {
      noOfGPS = p.second;
      usedCS = p.first;
    }
  }

  //if we didn't find any coordinate system we leave
  if(usedCS == CoordinateSystem())
    return NULL;

  //determine the dimensions of the bounding rect if the top left corner is exactly
  //the user choosen one and the bottom right corner is the one determined to include the user choosen one
  //given the determined cellsize
  const vector<RectCoord>& rcs = latLng2RC(asTlTrBrBl<vector<LatLngCoord> >(llrect), usedCS);
  auto rcpoly = Quadruple<RectCoord>(rcs);
  int noOfRows = int(std::ceil((rcpoly.tl.h - rcpoly.br.h)/minCellSize));
  int noOfCols = int(std::ceil((rcpoly.br.r - rcpoly.tl.r)/minCellSize));
  RectCoord tl(usedCS, min(rcpoly.tl.r, rcpoly.bl.r), max(rcpoly.tl.h, rcpoly.tr.h));
  RectCoord br(usedCS, tl.r + noOfCols*minCellSize, tl.h - noOfRows*minCellSize);

  map<string, GridPPtr> dsn2grid;
  set<string> uniqueDatasetNames;
  BOOST_FOREACH(CS2GMDS::value_type p, cs2gmds)
  {
    BOOST_FOREACH(const GridMetaData& gmd, p.second)
    {
      GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
      if(ci == gmd2gridProxies.end())
        continue;
      const GridProxies& agps = ci->second;
      BOOST_FOREACH(GridProxyPtr agp, agps)
      {
//        cout << "dsn: |" << agp->datasetName << "|" << endl;
        uniqueDatasetNames.insert(agp->datasetName);
      }
    }
  }

  //create all the grids needed
  BOOST_FOREACH(string dsn, uniqueDatasetNames)
  {
    dsn2grid[dsn] = GridPPtr(new GridP(dsn, noOfRows, noOfCols, minCellSize, tl.r, br.h, -9999, usedCS));
  }

  if(dsn2grid.empty())
    return NULL;

  string regionName;

  GridPPtr someGrid = dsn2grid.begin()->second;

  /*
  for(int r = 0, rs = someGrid->rows(); r < rs; r++)
  {
    for(int c = 0, cs = someGrid->cols(); c < cs; c++)
    {
      RectCoord cellCenter = someGrid->rcCoordAtCenter(r, c);
      double sgCellSize = someGrid->cellSize();

      LatLngCoord llcc = RC2latLng(cellCenter);

      BOOST_FOREACH(CS2GMDS::value_type p, cs2gmds)
      {
        CoordinateSystem cs = p.first;
        RectCoord rccc = latLng2RC(llcc, cs);

        BOOST_FOREACH(const GridMetaData& gmd, p.second)
        {
//          cout << "gmd: " << gmd.toCanonicalString() << " regionName: " << regionName << endl;
          regionName = gmd.regionName;
//          if(gmd.regionName == "uecker" ||
//             gmd.regionName == "sachsen" ||
//             gmd.regionName == "uelzen" ||
//             gmd.regionName == "brazil-sinop" ||
//             gmd.regionName == "brazil-campo-verde")
//            regionName = gmd.regionName;

          GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
          if(ci == gmd2gridProxies.end())
            continue;
          const GridProxies& agps = ci->second;
          if(agps.empty())
            continue;

          GridP* g = agps.front()->gridPtr();
          double gCellSize = g->cellSize();

          auto pos = g->rc2rowCol(rccc);

          //if any value lies farther outside than a cellwidth of the target grid, it can't intersect the closest
          //source cell thus we associate a no data value with the target grid value
          bool setNoDataValue =
              (pos.row < 0 && abs(pos.row*gCellSize) > sgCellSize) ||
              (pos.col < 0 && abs(pos.col*gCellSize) > sgCellSize) ||
              (pos.row > 0 && pos.isRowOutside && abs((pos.row - g->rows())*gCellSize > sgCellSize)) ||
              (pos.col > 0 && pos.isColOutside && abs((pos.col - g->cols())*gCellSize > sgCellSize));

          //if the position is in the cell outside, just associate the target grid value with the border value of
          //the source grid
          if(!setNoDataValue)
          {
            if(pos.row < 0)
              pos.row = 0;
            else if(pos.row > 0 && pos.isRowOutside)
              pos.row = g->rows() - 1;
            if(pos.col < 0)
              pos.col = 0;
            else if(pos.col > 0 && pos.isColOutside)
              pos.col = g->cols() - 1;
          }

          BOOST_FOREACH(GridProxyPtr agp, agps)
          {
            GridPPtr tg = dsn2grid[agp->datasetName];
            if(setNoDataValue)
              tg->setNoDataValueAt(r,c);
            else
              tg->setDataAt(r, c, agp->gridPtr()->dataAt(pos.row, pos.col));
          }
        }
      }
    }
  }
  */


  BOOST_FOREACH(CS2GMDS::value_type p, cs2gmds)
  {
    BOOST_FOREACH(const GridMetaData& gmd, p.second)
    {
      //          cout << "gmd: " << gmd.toCanonicalString() << " regionName: " << regionName << endl;
      regionName = gmd.regionName;
      //          if(gmd.regionName == "uecker" ||
      //             gmd.regionName == "sachsen" ||
      //             gmd.regionName == "uelzen" ||
      //             gmd.regionName == "brazil-sinop" ||
      //             gmd.regionName == "brazil-campo-verde")
      //            regionName = gmd.regionName;

      GMD2GPS::const_iterator ci = gmd2gridProxies.find(gmd);
      if(ci == gmd2gridProxies.end())
        continue;
      const GridProxies& agps = ci->second;
      if(agps.empty())
        continue;

      GridP* g = agps.front()->gridPtr();
      double gCellSize = g->cellSize();

      for(int r = 0, rows = someGrid->rows(); r < rows; r++)
      {
        for(int c = 0, cols = someGrid->cols(); c < cols; c++)
        {
          RectCoord cellCenter = someGrid->rcCoordAtCenter(r, c);
          double sgCellSize = someGrid->cellSize();

          LatLngCoord llcc = RC2latLng(cellCenter);

          CoordinateSystem cs = p.first;
          RectCoord rccc = latLng2RC(llcc, cs);

          auto pos = g->rc2rowCol(rccc);

          //if any value lies farther outside than a cellwidth of the target grid, it can't intersect the closest
          //source cell thus we associate a no data value with the target grid value
          bool setNoDataValue =
              (pos.row < 0 && abs(pos.row*gCellSize) > sgCellSize) ||
              (pos.col < 0 && abs(pos.col*gCellSize) > sgCellSize) ||
              (pos.row > 0 && pos.isRowOutside && abs((pos.row - g->rows())*gCellSize > sgCellSize)) ||
              (pos.col > 0 && pos.isColOutside && abs((pos.col - g->cols())*gCellSize > sgCellSize));

          //if the position is in the cell outside, just associate the target grid value with the border value of
          //the source grid
          if(!setNoDataValue)
          {
            if(pos.row < 0)
              pos.row = 0;
            else if(pos.row > 0 && pos.isRowOutside)
              pos.row = g->rows() - 1;
            if(pos.col < 0)
              pos.col = 0;
            else if(pos.col > 0 && pos.isColOutside)
              pos.col = g->cols() - 1;
          }

          BOOST_FOREACH(GridProxyPtr agp, agps)
          {
            GridPPtr tg = dsn2grid[agp->datasetName];
            if(setNoDataValue)
              tg->setNoDataValueAt(r,c);
            else
              tg->setDataAt(r, c, agp->gridPtr()->dataAt(pos.row, pos.col));
          }
        }
      }
    }
  }

  cout << "regionName: " << regionName << " noOfRows: " << noOfRows << " noOfCols: " << noOfCols << endl;

  auto vg = new VirtualGrid2(usedCS, RCRect(tl, br), minCellSize, noOfRows, noOfCols, dsn2grid);
  vg->setCustomId(regionName);

  return vg;
}

VirtualGrid2*
GridManager::virtualGridForGridMetaData(const GridMetaData& gmd,
                                        const Path& userSubPath)
{
	Path2GPS::const_iterator ci = _gmdMap.find(userSubPath);
  if(ci != _gmdMap.end())
  {
		GMD2GPS::const_iterator ci2 = ci->second.find(gmd);
    if(ci2 != ci->second.end())
    {
      NoVirtualGrid2* nvg =
          new NoVirtualGrid2(gmd.coordinateSystem,
														&(ci2->second),
														gmd.rcRect(), gmd.cellsize,
                            gmd.nrows, gmd.ncols);
			//cout << "gmd.regionName: " << gmd.regionName << endl;
			nvg->setName(gmd.regionName);
			nvg->setCustomId(gmd.regionName);
			return nvg;
		}
	}
{'fraction_non_alphanumeric': 0.16481415689658088, 'fraction_numerical': 0.008036001285760206, 'mean_word_length': 2.738772928526249, 'pattern_counts': {'":': 0, '<': 829, '<?xml version=': 0, '>': 366, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5409441', 'n_tokens_mistral': 21780, 'n_tokens_neox': 20131, 'n_words': 4669}