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

void ProcessPit_onepass(Array2D<elev_t> &dem, Array2D<label_t> &labels, std::queue<GridCellZ<elev_t> > &depressionQue, std::queue<GridCellZ<elev_t> > &traceQueue, GridCellZ_pq<elev_t> &priorityQueue, std::vector<std::map<label_t, elev_t> > &my_graph){
  while (!depressionQue.empty()){
    GridCellZ<elev_t> c = depressionQue.front();
    depressionQue.pop();

    for(int n=1;n<=8;n++){
      int nx = c.x+dx[n];
      int ny = c.y+dy[n];

      if(!dem.inGrid(nx,ny))
        continue;

      WatershedsMeet(labels(c.x,c.y),labels(nx,ny),dem(c.x,c.y),dem(nx,ny),my_graph);

      if(labels(nx,ny)!=0)
        continue;    

      labels(nx,ny) = labels(c.x,c.y);

      if (dem(nx,ny) > c.z) { //Slope cell
        traceQueue.emplace(nx,ny,dem(nx,ny));
      } else {                //Depression cell
        dem(nx,ny) = c.z;
        depressionQue.emplace(nx,ny,c.z);
      }
    }
  }
}

void ProcessTraceQue_onepass(Array2D<elev_t> &dem, Array2D<label_t> &labels, std::queue<GridCellZ<elev_t> > &traceQueue, GridCellZ_pq<elev_t> &priorityQueue, std::vector<std::map<label_t, elev_t> > &my_graph){
  while (!traceQueue.empty()){
    GridCellZ<elev_t> c = traceQueue.front();
    traceQueue.pop();

    bool bInPQ = false;
    for(int n=1;n<=8;n++){
      int nx = c.x+dx[n];
      int ny = c.y+dy[n];

      if(!dem.inGrid(nx,ny))
        continue;

      WatershedsMeet(labels(c.x,c.y),labels(nx,ny),dem(c.x,c.y),dem(nx,ny),my_graph);

      if(labels(nx,ny)!=0)
        continue;    
      
      //The neighbour is unprocessed and higher than the central cell
      if(c.z<dem(nx,ny)){
        traceQueue.emplace(nx,ny,dem(nx,ny));
        labels(nx,ny) = labels(c.x,c.y);
        continue;
      }

      //Decide  whether (nx, ny) is a true border cell
      if (!bInPQ) {
        bool isBoundary = true;
        for(int nn=1;nn<=8;nn++){
          int nnx = nx+dx[n];
          int nny = ny+dy[n];
          if(!dem.inGrid(nnx,nny))
            continue;
          if (labels(nnx,nny)!=0 && dem(nnx,nny)<dem(nx,ny)){
            isBoundary = false;
            break;
          }
        }
        if(isBoundary){
          priorityQueue.push(c);
          bInPQ = true;
        }
      }
    }
  }
}

static inline TA_Setup_Vars TerrainSetup(const Array2D<T> &elevations, const int x, const int y, const float zscale){
  TA_Setup_Vars tsv;
  tsv.a=tsv.b=tsv.c=tsv.d=tsv.e=tsv.f=tsv.g=tsv.h=tsv.i=elevations(x,y);
  if(elevations.inGrid(x-1,y-1) && elevations(x-1,y-1)!=elevations.noData()) tsv.a = elevations(x-1,y-1);
  if(elevations.inGrid(x-1,y  ) && elevations(x-1,y  )!=elevations.noData()) tsv.d = elevations(x-1,y  );
  if(elevations.inGrid(x-1,y+1) && elevations(x-1,y+1)!=elevations.noData()) tsv.g = elevations(x-1,y+1);
  if(elevations.inGrid(x  ,y-1) && elevations(x,  y-1)!=elevations.noData()) tsv.b = elevations(x,  y-1);
  if(elevations.inGrid(x  ,y+1) && elevations(x,  y+1)!=elevations.noData()) tsv.h = elevations(x,  y+1);
  if(elevations.inGrid(x+1,y-1) && elevations(x+1,y-1)!=elevations.noData()) tsv.c = elevations(x+1,y-1);
  if(elevations.inGrid(x+1,y  ) && elevations(x+1,y  )!=elevations.noData()) tsv.f = elevations(x+1,y  );
  if(elevations.inGrid(x+1,y+1) && elevations(x+1,y+1)!=elevations.noData()) tsv.i = elevations(x+1,y+1);

  tsv.a *= zscale;
  tsv.b *= zscale;
  tsv.c *= zscale;
  tsv.d *= zscale;
  tsv.e *= zscale;
  tsv.f *= zscale;
  tsv.g *= zscale;
  tsv.h *= zscale;
  tsv.i *= zscale;

  return tsv;
}

label_t GetNewLabelZhou(
  int x,
  int y,
  label_t &current_label,
  uint8_t edge,
  const Array2D<elev_t> &dem,
  const Array2D<label_t> &labels
){
  if(labels(x,y)!=0)
    return labels(x,y);

  for(int n=1;n<=8;n++){
    int nx = x+dx[n];
    int ny = y+dy[n];
    if(!dem.inGrid(nx,ny))
      continue;
    if(labels(nx,ny)!=0 && dem(nx,ny)<=dem(x,y))
      return labels(nx,ny);
  }

  return current_label++;
}

void Zhou2015Labels(
  Array2D<elev_t>                         &dem,
  Array2D<label_t>                        &labels,
  std::vector<std::map<label_t, elev_t> > &my_graph,
  uint8_t edge,
  bool    flipH,
  bool    flipV
){
  std::queue<GridCellZ<elev_t> > traceQueue;
  std::queue<GridCellZ<elev_t> > depressionQue;

  label_t current_label = 2;

  labels.setAll(0);

  GridCellZ_pq<elev_t> priorityQueue;

  for(int32_t x=0;x<dem.width();x++){
    const int height = dem.height()-1;
    priorityQueue.emplace(x,0,     dem(x,0     ));
    priorityQueue.emplace(x,height,dem(x,height));
  }

  for(int32_t y=1;y<dem.height()-1;y++){
    const int width = dem.width()-1;
    priorityQueue.emplace(0,    y,dem(0,    y));
    priorityQueue.emplace(width,y,dem(width,y));
  }

  while (!priorityQueue.empty()){
    GridCellZ<elev_t> c = priorityQueue.top();
    priorityQueue.pop();

    auto my_label = labels(c.x,c.y) = GetNewLabelZhou(c.x,c.y,current_label,edge,dem,labels);

    for(int n=1;n<=8;n++){
      int nx = c.x+dx[n];
      int ny = c.y+dy[n];

      if (!dem.inGrid(nx,ny))
        continue;

      WatershedsMeet(my_label,labels(nx,ny),dem(c.x,c.y),dem(nx,ny),my_graph);

      if(labels(nx,ny)!=0)
        continue;

      labels(nx,ny) = labels(c.x,c.y);

      if(dem(nx,ny)<=c.z){ //Depression cell
        dem(nx,ny) = c.z;
        depressionQue.emplace(nx,ny,c.z);
        ProcessPit_onepass(dem,labels,depressionQue,traceQueue,priorityQueue,my_graph);
      } else {          //Slope cell
        traceQueue.emplace(nx,ny,dem(nx,ny));
      }     
      ProcessTraceQue_onepass(dem,labels,traceQueue,priorityQueue,my_graph);
    }
  }

  //Connect the DEM's outside edges to Special Watershed 1. This requires
  //knowing whether the tile has been flipped toe snure that we connect the
  //correct edges.
  if( ((edge & GRID_TOP)    && !flipV) || ((edge & GRID_BOTTOM) && flipV) )
    for(int32_t x=0;x<labels.width();x++)
      WatershedsMeet(labels(x,0),(label_t)1,dem(x,0),dem(x,0),my_graph);

  if( ((edge & GRID_BOTTOM) && !flipV) || ((edge & GRID_TOP)    && flipV) ){
    int bottom_row = labels.height()-1;
    for(int32_t x=0;x<labels.width();x++)
      WatershedsMeet(labels(x,bottom_row),(label_t)1,dem(x,bottom_row),dem(x,bottom_row),my_graph);
  }

  if( ((edge & GRID_LEFT)  && !flipH) || ((edge & GRID_RIGHT) && flipH) )
    for(int32_t y=0;y<labels.height();y++)
      WatershedsMeet(labels(0,y),(label_t)1,dem(0,y),dem(0,y),my_graph);  

  if( ((edge & GRID_RIGHT) && !flipH) || ((edge & GRID_LEFT)  && flipH) ){
    int right_col = labels.width()-1;
    for(int32_t y=0;y<labels.height();y++)
      WatershedsMeet(labels(right_col,y),(label_t)1,dem(right_col,y),dem(right_col,y),my_graph);
  }

  my_graph.resize(current_label);
}

void dinf_upslope_area(
  const Array2D<T> &flowdirs,
  Array2D<U> &area
){
  Array2D<int8_t> dependency;
  std::queue<GridCell> sources;
  ProgressBar progress;

  std::cerr<<"\nA D-infinity Upslope Area"<<std::endl;
  std::cerr<<"C Tarboton, D.G. 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research. Vol. 33. pp 309-319."<<std::endl;

  std::cerr<<"p Setting up the dependency matrix..."<<std::endl;
  dependency.resize(flowdirs);
  dependency.setAll(0);

  std::cerr<<"p Setting up the area matrix..."<<std::endl;
  area.resize(flowdirs);
  area.setAll(0);
  area.setNoData(dinf_NO_DATA);

  bool has_cells_without_flow_directions=false;
  std::cerr<<"p Calculating dependency matrix & setting noData() cells..."<<std::endl;
  progress.start( flowdirs.size() );

  ///////////////////////
  //Calculate the number of "dependencies" each cell has. That is, count the
  //number of cells which flow into each cell.

  #pragma omp parallel for reduction(|:has_cells_without_flow_directions)
  for(int y=0;y<flowdirs.height();y++){
    progress.update( y*flowdirs.width() );
    for(int x=0;x<flowdirs.width();x++){
      //If the flow direction of the cell is NoData, mark its area as NoData
      if(flowdirs.isNoData(x,y)){
        area(x,y)       = area.noData();
        dependency(x,y) = 9;  //TODO: This is an unnecessary safety precaution. This prevents the cell from ever being enqueued (an unnecessary safe guard? TODO)
        continue;             //Only necessary if there are bugs below (TODO)
      }

      //If the cell has no flow direction, note that so we can warn the user
      if(flowdirs(x,y)==NO_FLOW){
        has_cells_without_flow_directions=true;
        continue;
      }

      //TODO: More explanation of what's going on here
      int n_high, n_low;
      int nhx,nhy,nlx,nly;
      where_do_i_flow(flowdirs(x,y),n_high,n_low);
      nhx=x+dinf_dx[n_high];
      nhy=y+dinf_dy[n_high];
      if(n_low!=-1){
        nlx = x+dinf_dx[n_low];
        nly = y+dinf_dy[n_low];
      }
      if( n_low!=-1 && flowdirs.inGrid(nlx,nly) && flowdirs(nlx,nly)!=flowdirs.noData() )
        dependency(nlx,nly)++;
      if( flowdirs.inGrid(nhx,nhy) && flowdirs(nhx,nhy)!=flowdirs.noData() )
        dependency(nhx,nhy)++;
    }
  }
  std::cerr<<"t Succeeded in = "<<progress.stop()<<" s"<<std::endl;

  if(has_cells_without_flow_directions)
    std::cerr<<"W \033[91mNot all cells had defined flow directions! This implies that there will be digital dams!\033[39m"<<std::endl;

  ///////////////////////
  //Find those cells which have no dependencies. These are the places to start
  //the flow accumulation calculation.

  std::cerr<<"p Locating source cells..."<<std::endl;
  progress.start( flowdirs.size() );
  for(int y=0;y<flowdirs.height();y++){
    progress.update( y*flowdirs.width() );
    for(int x=0;x<flowdirs.width();x++)
      if(flowdirs(x,y)==flowdirs.noData())
        continue;
      else if(flowdirs(x,y)==NO_FLOW)
        continue;
      else if(dependency(x,y)==0)
        sources.emplace(x,y);
  }
  std::cerr<<"t Source cells located in = "<<progress.stop()<<" s"<<std::endl;





  ///////////////////////
  //Calculate the flow accumulation by "pouring" a cell's flow accumulation
  //value into the cells below it, as indicated by the D-infinite flow routing
  //method.

  std::cerr<<"p Calculating up-slope areas..."<<std::endl;
  progress.start( flowdirs.numDataCells() );
  long int ccount=0;
  while(sources.size()>0){
    auto c = sources.front();
    sources.pop();

//.........这里部分代码省略.........

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();
}