C++ Array2D::isEdgeCell方法代码示例

static int d8_FlowDir(const Array2D<T> &elevations, const int x, const int y){
  T minimum_elevation = elevations(x,y);
  int flowdir         = NO_FLOW;

  if (elevations.isEdgeCell(x,y)){
    if(elevations.isTopLeft(x,y))
      return 2;
    else if(elevations.isBottomLeft(x,y))
      return 8;
    else if(elevations.isTopRight(x,y))
      return 4;
    else if(elevations.isBottomRight(x,y))
      return 6;
    else if(elevations.isLeftCol(x,y))
      return 1;
    else if(elevations.isRightCol(x,y))
      return 5;
    else if(elevations.isTopRow(x,y))
      return 3;
    else if(elevations.isBottomRow(x,y))
      return 7;
  }

  /*NOTE: Since the very edges of the DEM are defined to always flow outwards,
  if they have defined elevations, it is not necessary to check if a neighbour
  is IN_GRID in the following
  NOTE: It is assumed that the no_data datum is an extremely negative
  number, such that all water which makes it to the edge of the DEM's region
  of defined elevations is sucked directly off the grid, rather than piling up
  on the edges.*/
  for(int n=1;n<=8;n++)
    if(
      elevations(x+dx[n],y+dy[n])<minimum_elevation
      || (elevations(x+dx[n],y+dy[n])==minimum_elevation
            && flowdir>0 && flowdir%2==0 && n%2==1) //TODO: What is this modulus stuff for?
    ){
      minimum_elevation=elevations(x+dx[n],y+dy[n]);
      flowdir=n;
    }

  return flowdir;
}

void FindFlats(
  const Array2D<T>   &elevations,
  Array2D<int8_t>    &flats
){
  flats.resize(elevations);
  flats.setNoData(FLAT_NO_DATA);

  ProgressBar progress;

  progress.start( elevations.size() );

  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++)
  for(int x=0;x<elevations.width();x++){
    if(elevations.isNoData(x,y)){
      flats(x,y) = FLAT_NO_DATA;
      continue;
    }

    if(elevations.isEdgeCell(x,y)){
      flats(x,y) = NOT_A_FLAT;
      continue;
    }

    //We'll now assume that the cell is a flat unless proven otherwise
    flats(x,y) = IS_A_FLAT;

    for(int n=1;n<=8;n++){
      const int nx = x+dx[n];
      const int ny = y+dy[n];
      if(elevations(nx,ny)<elevations(x,y) || elevations.isNoData(nx,ny)){
        flats(x,y) = NOT_A_FLAT;
        break;
      }
    }

    //We handled the base case just above the for loop
  }

  RDLOG_TIME_USE<<"Succeeded in = "<<progress.stop()<<" s";
}

static float dinf_FlowDir(const Array2D<T> &elevations, const int x, const int y){
  //Ensure that flow is pulled off the edge of the grid
  if (elevations.isEdgeCell(x,y)){
    if(x==0 && y==0)
      return 3*M_PI/4;  //D8: 2
    else if(x==0 && y==elevations.height()-1)
      return 5*M_PI/4;  //D8: 8
    else if(x==elevations.width()-1 && y==0)
      return 1*M_PI/4;  //D8: 4
    else if(x==elevations.width()-1 && y==elevations.height()-1)
      return 7*M_PI/4;  //D8: 6
    else if(x==0)
      return 4*M_PI/4;  //D8: 1
    else if(x==elevations.width()-1)
      return 0*M_PI/4;  //D8: 5
    else if(y==0)
      return 2*M_PI/4;  //D8: 3
    else if(y==elevations.height()-1)
      return 6*M_PI/4;  //D8: 7
  }

  int    nmax = -1;
  double smax = 0;
  double rmax = 0;

  //I am not on the edge of the grid. All my neighbours can be examined.

  for(int n=0;n<8;n++){
    //Is is assumed that cells with a value of NoData have very negative
    //elevations with the result that they draw flow off of the grid.

    //Choose elevations based on Table 1 of Tarboton (1997), Barnes TODO
    const double e0 = elevations(x,y);
    const double e1 = elevations(x+dx_e1[n],y+dy_e1[n]);
    const double e2 = elevations(x+dx_e2[n],y+dy_e2[n]);

    //TODO: Assumes that the width and height of grid cells are equal and scaled
    //to 1.
    const double d1 = 1;
    const double d2 = 1;

    const double s1 = (e0-e1)/d1;
    const double s2 = (e1-e2)/d2;
    double r        = atan2(s2,s1);

    double s;

    if(r<0){
      r = 0;
      s = s1;
    } else if(r>atan2(d2,d1)){
      r = atan2(d2,d1); //TODO: This is a constant
      s = (e0-e2)/sqrt(d1*d1+d2*d2);
    } else {
      s = sqrt(s1*s1+s2*s2);
    }

    if(s>smax){
      smax = s;
      nmax = n;
      rmax = r;
    }
  }

  double rg = NO_FLOW;
  if(nmax!=-1)
    rg = (af[nmax]*rmax+ac[nmax]*M_PI/2);

  return rg;
}

void FM_Freeman(
  const Array2D<E> &elevations,
  Array3D<float> &props,
  const double xparam
){
  RDLOG_ALG_NAME<<"Freeman (1991) Flow Accumulation (aka MFD, MD8)";
  RDLOG_CITATION<<"Freeman, T.G., 1991. Calculating catchment area with divergent flow based on a regular grid. Computers & Geosciences 17, 413–422.";
  RDLOG_CONFIG<<"p = "<<xparam;

  props.setAll(NO_FLOW_GEN);
  props.setNoData(NO_DATA_GEN);

  ProgressBar progress;
  progress.start(elevations.size());

  #pragma omp parallel for collapse(2)
  for(int y=0;y<elevations.height();y++)
  for(int x=0;x<elevations.width();x++){
    ++progress;

    if(elevations.isNoData(x,y)){
      props(x,y,0) = NO_DATA_GEN;
      continue;
    }

    if(elevations.isEdgeCell(x,y))
      continue;

    const E e    = elevations(x,y);

    double C = 0;
    for(int n=1;n<=8;n++){
      const int nx = x+dx[n];
      const int ny = y+dy[n];

      if(!elevations.inGrid(nx,ny))
        continue;
      if(elevations.isNoData(nx,ny)) //TODO: Don't I want water to drain this way?
        continue;

      const E ne = elevations(nx,ny);

      if(ne<e){
        const double rise = e-ne;
        const double run  = dr[n];
        const double grad = rise/run;
        const auto cval   = std::pow(grad,xparam);
        props(x,y,n)      = cval;
        C                += cval;
      }
    }

    if(C>0){
      props(x,y,0) = HAS_FLOW_GEN;

      C = 1/C; //TODO

      for(int n=1;n<=8;n++){
        auto &this_por = props(x,y,n);
        if(this_por>0)
          this_por *= C;
        else
          this_por = 0;
      }
    }
  }
  progress.stop();
}