grid.cpp

Go to the documentation of this file.

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

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <algorithm> // for sorting
#include <sstream>
#include <valarray>
#include <time.h> // for seed of random number generator
#include <math.h>
#include "r250.h"
#include "grid.h"
//---------------------------------------------------------------------------

float RandG(int, int) {return float(dr250());}

bool operator<(const xyC& a, const xyC& b) // new 20051222
{
  return a.height < b.height;
};

std::string itoa(int value) {
  
  enum { kMaxDigits = 35 };
  
  std::string buf;
  
  buf.reserve( kMaxDigits ); // Pre-allocate enough space.
  
  int quotient = value;
  
  // Translating number to string with base:  
  do {  
        buf += "0123456789"[ std::abs(quotient%10) ];   
        quotient /= 10; 
  } while ( quotient );
  
  // Append the negative sign for base 10  
  if ( value < 0) buf += '-';
  
  std::reverse( buf.begin(), buf.end() );        
  return buf;        
}


TGrid::TGrid (void)
{ 
  time_t currentTime;
  int seed = time(&currentTime); // seed for random number generator, should be made variable

  r250_init(29); // start of the random number generator
  
  clearArrays(); // set all arrays to zero in a new grid
}//konstruktor
//---------------------------------------------------------------------------
TGrid::~TGrid ()
{

}//Destruktor
//---------------------------------------------------------------------------

void TGrid::clearArrays(void) // new 20051222
{
  for (int x=0; x<=N; x++)
    for (int y=0; y<=N; y++)
    {
       CellArray[x][y]=0;
       CellArrayF[x][y]=0;
       CellArrayL[x][y]=0;
       CellCount[x][y].Position=0;
           CellCount[x][y].NCellsAll=0;
    }//for y
}
//---------------------------------------------------------------------------
inline float TGrid::f3 (float delta, float x0, float x1, float x2)
{
  return(x0+x1+x2)/3.+delta*RandG (0,1);                       //RandG (0,1);
} 
//---------------------------------------------------------------------------
inline float TGrid::f4 (float delta, float x0, float x1, float x2, float x3)
{
  return(x0+x1+x2+x3)/4.+delta*RandG (0,1);
}
//---------------------------------------------------------------------------
void TGrid::createMountain (int Iterations,int Sigma,float Hurst,bool Addition)
{
// set the initial random corners
float delta=Sigma;
CellArrayF[0][0]=delta*RandG (0,1);
CellArrayF[0][N]=delta*RandG (0,1);
CellArrayF[N][0]=delta*RandG (0,1);
CellArrayF[N][N]=delta*RandG (0,1);
int D=N;
int d=N/2;

for (int stage=1; stage<=Iterations; stage++)
{
  delta=delta*pow(0.5,0.5*Hurst); //going from Type I to type II
                                                                  //interpolate and offset points
  for (int x=d; x<=(N-d);x+=D)
        for (int y=d; y<=(N-d);y+=D)
          CellArrayF[x][y]=f4(delta,CellArrayF[x+d][y+d],CellArrayF[x+d][y-d],CellArrayF[x-d][y+d],CellArrayF[x-d][y-d]);

          //displace other points also, if needed
  if (Addition==true)
  {
        for (int x=0; x<=N; x+=D)
          for (int y=0; y<=N; y+=D)
                CellArrayF[x][y]+= delta*RandG (0,1);
  }//if addition
  
  //going from grid type II to type I
  delta *= pow(0.5,0.5*Hurst);
  
  //interpolate and offset boundary grid points
  for (int x=d; x<=(N-d); x+=D)
  {
        CellArrayF[x][0]=f3(delta,CellArrayF[x+d][0],CellArrayF[x-d][0],CellArrayF[x][d]);
        CellArrayF[x][N]=f3(delta,CellArrayF[x+d][N],CellArrayF[x-d][N],CellArrayF[x][N-d]);
        CellArrayF[0][x]=f3(delta,CellArrayF[0][x+d],CellArrayF[0][x-d],CellArrayF[d][x]);
        CellArrayF[N][x]=f3(delta,CellArrayF[N][x+d],CellArrayF[N][x-d],CellArrayF[N-d][x]);
  }//for
  
   //interpolate and offset interior Points
  for (int x=d; x<=(N-d);x+=D)
        for (int y=D; y<=(N-d);y+=D)
          CellArrayF[x][y]=f4(delta,CellArrayF[x][y+d],CellArrayF[x][y-d],CellArrayF[x+d][y],CellArrayF[x-d][y]);

  for (int x=D; x<=(N-d);x+=D)
        for (int y=d; y<=(N-d);y+=D)
          CellArrayF[x][y]=f4(delta,CellArrayF[x][y+d],CellArrayF[x][y-d],CellArrayF[x+d][y],CellArrayF[x-d][y]);
  
  //displace other points if needed
  if (Addition==true)
  {
        for (int x=0; x<=N; x+=D)
          for (int y=0; y<=N; y+=D)
                CellArrayF[x][y]+= delta*RandG (0,1);

        for (int x=d; x<=(N-d); x+=D)
          for (int y=d; y<=(N-d); y+=D)
                CellArrayF[x][y]+= delta*RandG (0,1);
  }//if Addition
  D=D/2;
  d=d/2;
}//for stage
} // end createMountain

//---------------------------------------------------------------------------
void TGrid::createWadi (float Grad, bool Parabel)
{
  int Center=N/2;
  int DiffY; // distance from centre of wadi cross-section
  float gradient = Grad?tan(Grad/360.*2*3.14159):0.0; // transformation of degrees to rad and to a slope angle defined in terms of the length of the adjacent leg (Ankathete) ; new 20061016
  
  if (Parabel==false)
  {
        for (int x=0;x<=N;x++)
          for (int y=0;y<=N;y++)
          {
                DiffY = Center-y;  // distance from centre of wadi cross-section

                /*== create a V-shaped curve with adjustable minimum ==*/
                CellArrayL [x][y]= (-gradient*26 + abs(DiffY)*gradient)*5; // setting elevation (m)             
                if (abs(DiffY)<=25) // cut the V-curve at zero
                { CellArrayL [x][y]=0.0;
                  CellCount [x][y].Position='W';// set Wadi habitat
                }
                if (Grad)
                {               
                  if (CellArrayL [x][y]>=100)   
                  {   CellArrayL [x][y] =100;   // cut the V-curve at 100 m
                        CellCount [x][y].Position='P';//set Plateau habitat
                  }
                  if (DiffY>25 && CellArrayL [x][y]<100) CellCount [x][y].Position='N';//set North facing habitat
                  if (DiffY<-25 && CellArrayL [x][y]<100) CellCount [x][y].Position='S';//set South facing habitat
                }
                else // zero degrees, flat landscape; corrected 20061031
                {
                  if (DiffY>133) CellCount [x][y].Position='P'; // set plateau habitat
                  else if (DiffY>25) CellCount [x][y].Position='N';// set North facing slope habitat
                  else if (DiffY>=-25) CellCount [x][y].Position='W'; // set Wadi habitat
                  else if (DiffY>-134) CellCount [x][y].Position='S';//set South facing slope habitat
                  else CellCount [x][y].Position='P'; // set plateau habitat
                }
          }//for y

        /*== smoothe the edges of the wadiscape parallel to cross-section with a moving average ==*/
        float* MemArrayL = new float [(N+1)*(N+1)]; // temporary grid for smoothed elevation; new 20051222
        int smoother = 10; // length of moving average window
        float counter = 0.0; // changed 20060111

        for (int x=0;x<=N;x++) // go parallel to cross-section
          for (int y=smoother;y<=N-smoother;y++) // changed 20060110
          {
                // add elevation in moving average window
                counter=0;
                for (int iz=-smoother; iz<=smoother; iz++)
                  counter+= CellArrayL[x][y+iz];
                
                // put average into temporary grid
        MemArrayL[x*(N+1)+y] = (1.0/(2.0*smoother+1.0))*counter;          
          }  //for x, y
                
        for (int x=0;x<=N;x++)
          for (int y=smoother;y<=N-smoother;y++) // changed 20060110
                CellArrayL[x][y]=MemArrayL[x*(N+1)+y]; // copy back the values
        
        delete(MemArrayL); // new 20051222
  }//if not parabola
  
  /*== create a wadi with parabolic cross section ==*/
  else if (Parabel==true)
  {
    for (int x=0;x<=N;x++)
      for (int y=0;y<=N;y++)
      {
                DiffY = Center == y ? 1 : Center - y;
                CellArrayL [x][y]=(gradient/-100.0)*(DiffY*DiffY);
      }//for x,y
  }// else if parabola
        
}//createWadi


//---------------------------------------------------------------------------
void TGrid::calculateElevation ()
{
  /*== add the elevation of all cells ==*/
  float sum_elevation = 0.0;
  for(int i = N/2-150; i<N/2+150; i++)
        for(int j = N/2-150; j<N/2+150; j++)
          sum_elevation += CellArrayF [i][j];
          
  /*== subtract the mean elevation from each cell's elevation ==*/
          for(int i = 0; i<=N; i++)
                for(int j = 0; j<=N; j++)
                  CellArrayF [i][j] -= sum_elevation/float(300*300); // changed 20070215

  /*== drape fractal landscape over wadi landscape ==*/
  for (int i=0;i<=N;i++)
    for (int j=0;j<=N;j++)
          CellArray [i][j] = CellArrayL [i][j] + CellArrayF [i][j]; // The relation of fractal to wadi elevation could be changed here.   
  
  /*== smoothe wadiscape in 2 dimensions; new 20060110 ==*/
  float* MemArray = new float [(N+1)*(N+1)]; // temporary grid for smoothed elevation; new 20051222 
  int smoother = 2; // length of moving average window
  float counter; // changed 20060111
   
  for (int x=0; x<=N;x++)  // changed 20060110
    for (int y=0; y<=N;y++)  // changed 20060110
    {
          int a=0;
          // add elevation in moving average window
          counter=0; // reset
          for (int iz = -smoother; iz<= smoother; iz++) 
                for (int is = -smoother; is<= smoother; is++)
          counter += CellArray[x+is][y+iz];
                  
          // put average into temporary grid
          MemArray[x*(N+1)+y]=(1.0/((2.0*smoother+1.0)*(2.0*smoother+1.0)))*counter; // new 20051222
  }  //for x, y
  
  for (int x=smoother;x<=N-smoother;x++) // changed 20060110
    for (int y=smoother;y<=N-smoother;y++) // changed 20060110
          CellArray[x][y]=MemArray[x*(N+1)+y];  // copy back the values
  
  delete MemArray; // new 20060110 
  
}//calculateElevation ()

//---------------------------------------------------------------------------
// neu 20051222
float flowinfo::getTotalContributing (flowinfo* receivingCell)
{ 
  if (checked == false)
  {
        for (int c=0; c<contributing; c++) // from all contributing cells
          if (receivingCell != nbsFrom[c])
                total += nbsFrom[c]->getTotalContributing(this); // get their number of contributing cells
        checked = true;
  }
  return total*share;
}

//---------------------------------------------------------------------------
// neu 20051206
void flowinfo::setInflow (flowinfo* pCell)
{  nbsFrom[contributing++] = pCell;
}

//---------------------------------------------------------------------------
void TGrid::countHigherCellsGIS (void) // neu 20051222
{
  /*== max. variation to add to neighbour cells with the same height ==*/
  // find out the lowest and highest point
  float maxHeight = CellArray[N/2][N/2]; // new 20051222
  float minHeight = CellArray[N/2][N/2]; // new 20051222
  for (int x=0; x<=N; x++)
        for (int y=0; y<=N; y++)
        {       
          maxHeight = CellArray[x][y]>maxHeight?CellArray[x][y]:maxHeight;
          minHeight = CellArray[x][y]<minHeight?CellArray[x][y]:minHeight;
        }
  
  const float smalldiff = (maxHeight-minHeight)/100.0; // changed 20051222
  

  /*== determine flow accumulation ==*/
  for (int x=0; x<=N; x++)
        for (int y=0; y<=N; y++)
        {       
          float a[8]; // angle from top of centre cell to top of each nb cell; water is assumed to flow down to the cell with the greatest negative angle.
          int nnbs = 0; // ID of neighbours
          flowinfo* myCell[8]; // addresses of neighbours
          
          for (int i=-1; i<=1; i++)
                for (int j=-1; j<=1; j++)
                { if (x+i < 0 || x+i > N) // avoid edge
                        continue;
                  if ( y+j < 0 || y+j > N) // avoid edge
                        continue;
                  if (!i && !j) // exempt self
                        continue;
                  
                  /* if direct neighbour cells have the same height, add a little variation to avoid endless loops in the recursive function getTotalContributing. */
//                if (CellArray[x+i][y+j] == CellArray[x][y])
//                      CellArray[x][y] += (RandG(0,1)-0.5)*2.0*smalldiff;
                  
                  // calculate angle in triangle
                  float ak = CellArray[x+i][y+j]-CellArray[x][y]; // adjacent leg (Ankathete)
                  float gk = sqrt(i*i + j*j); // opposite leg (Gegenkathete)
                  if (ak) a[nnbs] = atan(gk/ak); // remember all angles
                  else a[nnbs] = -3.14159/2.0; // same elevation, set a[] to zero; this case should never occur because I add a smalldiff when two neighbouring cells have the same elevation (see above)
                  myCell[nnbs] = &dir[x+i][y+j]; // store addresses of info of each nb

                  nnbs++; // count neighbours - important at edges
                } // end for i,j
                  
                float aMin = 0; // -90∞
                int ties=-1; 
                flowinfo* nbsTo[4]; // store info about nbs, there can be no more than 4 nbs with the same angle
                for (int n=0; n<nnbs; n++) // look for smallest angle downwards
                { if (a[n]<aMin) // if downward and greater than current steepest angle; changed 20070215
                  { aMin=a[n];  
                        ties=0; // if a steeper angle was found, start over at index zero
                        nbsTo[ties] = myCell[n]; // nbsTo contains the addresses of cells where water flows to
                  } // split flow if two or more nb cells have the same distance and height
                  else if (a[n] == aMin)
                  { ties++;
                        nbsTo[ties] = myCell[n]; // store in tiesMem the addresses of cells where water flows to
                  }
                } // end for n
                
                
                for (int t=0; t<ties+1; t++)                    
                  if (nbsTo[t] != &dir[x][y]) // exclude self
                   nbsTo[t]->setInflow(&dir[x][y]); // tell the lowest nbs where the water comes from
                  
                if(ties>-1) dir[x][y].setShare(1.0/float(ties+1)); // if there were ties, split the water
  } //end for x,y
  
  /* *****************************************
        There were segmentation errors (SIGSEGV) when the grid was checked in the usual row-by-column way. Because getTotalContributing() is a recursive function, it could hit a very low cell and end up collecting data from 10'000+ cells. Now the whole grid is sorted by height and the contributing cells are checked from the highest to the lowest.
        ***************************************** */
  xyC* pDir = new xyC[(N+1)*(N+1)];
  for (int x=0; x<=N; x++)
        for (int y=0; y<=N; y++)
        {
          pDir[x*(N+1)+y].col = x;
          pDir[x*(N+1)+y].row = y;
          pDir[x*(N+1)+y].height = CellArray[x][y];
        }
          
        std::sort(pDir, pDir+(N+1)*(N+1));
  /******************************************/
 
  for (long w=(N+1)*(N+1)-1; w>=0; w--)
  { 
        int x = pDir[w].col; int y = pDir[w].row;
        float TC = dir[x][y].getTotalContributing(&dir[x][y]);
        CellCount [x][y].NCellsAll = TC<1.0?1.0:TC; // peaks with shared flow should be shown as one patch
  }
  delete(pDir);
  
}


//---------------------------------------------------------------------------
void TGrid::outputTestVariable(void)
{
  char* fn = "TestVariable.txt";
  std::ofstream Ausgabe;
  Ausgabe.open(fn);
  for (int i=0;i<=N; i++)
  {     for (int j=0; j<=N; j++)
        {
          Ausgabe<< CellArrayF[i][j];
          if (j<N)  Ausgabe  << "\t";
        }
        if (i<N)  Ausgabe << std::endl;
  }
  
  Ausgabe.close();
}

//---------------------------------------------------------------------------
void TGrid::outputElevation(void)
{
  char* fn = "Elevation.txt";
  std::ofstream Ausgabe;
  Ausgabe.open(fn);
  for (int i=0;i<=N; i++)
  {     for (int j=0; j<=N; j++)
  {
          Ausgabe<< CellArray[i][j];
        if (j<N)  Ausgabe  << "\t";
  }
        if (i<N)  Ausgabe << std::endl;
  }

  Ausgabe.close();
}

//---------------------------------------------------------------------------
void TGrid::outputContributing(void)
{
  char* fn = "Contributing.txt";
  std::ofstream Ausgabe;
  Ausgabe.open(fn);
  for (int i=0;i<=N; i++)
  {     for (int j=0; j<=N; j++)
  {
        Ausgabe << (CellCount[i][j]).NCellsAll;
        if (j<N)  Ausgabe  << "\t";
  }
        if (i<N)  Ausgabe << std::endl;
  }
  Ausgabe.close();
}

//---------------------------------------------------------------------------
void TGrid::writeFile(float slope, int r)
{
  std::string filename = "Wadiscape_slope_";
  std::string tag = ".txt";
  filename += itoa(int(slope)) ;
  filename += "_";
  filename += itoa(r);
  filename += tag;
  
  int center=N/2;
  
  std::ofstream DataFile;
  DataFile.open(filename.c_str()) ;
  if (!DataFile)
  {
        //ˆffnen ist gescheitert =>Fehlerbehandlung
        std::cerr<<"Die Ausgabedatei konnte nicht erfolgreich"
        " geˆffnet werden."<< std::endl;
        exit(0);
  }//if kein DataFile
  int xx = 150;
  for (int i=center-xx;i<center+xx;i++)
  {
    for (int j=center-xx;j<center+xx;j++)
    {
      DataFile<<i-(center-xx)<< "\t" << j-(center-xx)<< "\t"<< CellArray[i][j]<<"\t"; // X and Y coordinates
      DataFile<<CellCount[i][j].Position<<"\t"; // habitat
          DataFile<<CellCount[i][j].NCellsAll<<"\t"; //
                DataFile<<slope; // slope angle in degrees
      DataFile<<"\n";
    }//for j
  }//for i
  DataFile.close ();

}

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

---