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

void TA_CTI(
  const Array2D<T> &flow_accumulation,
  const Array2D<U> &riserun_slope,
        Array2D<V> &result
){
  Timer timer;

  RDLOG_ALG_NAME<<"d8_CTI";

  if(flow_accumulation.width()!=riserun_slope.width() || flow_accumulation.height()!=riserun_slope.height())
    throw std::runtime_error("Couldn't calculate CTI! The input matricies were of unequal dimensions!");

  RDLOG_PROGRESS<<"Setting up the CTI matrix..."<<std::flush;
  result.resize(flow_accumulation);
  result.setNoData(-1);  //Log(x) can't take this value of real inputs, so we're good
  RDLOG_PROGRESS<<"succeeded.";

  RDLOG_PROGRESS<<"Calculating CTI..."<<std::flush;
  timer.start();
  #pragma omp parallel for collapse(2)
  for(int x=0;x<flow_accumulation.width();x++)
    for(int y=0;y<flow_accumulation.height();y++)
      if(flow_accumulation(x,y)==flow_accumulation.noData() || riserun_slope(x,y)==riserun_slope.noData())
        result(x,y)=result.noData();
      else
        result(x,y)=log( (flow_accumulation(x,y)/flow_accumulation.getCellArea()) / (riserun_slope(x,y)+0.001) );
  RDLOG_TIME_USE<<"succeeded in "<<timer.stop()<<"s.";
}

void d8_flow_directions(
  const Array2D<T> &elevations,
        Array2D<U> &flowdirs
){
  ProgressBar progress;

  std::cerr<<"A D8 Flow Directions"<<std::endl;
  std::cerr<<"C TODO"<<std::endl;

  std::cerr<<"p Setting up the flow directions matrix..."<<std::endl;
  flowdirs.resize(elevations);
  flowdirs.setAll(NO_FLOW);
  flowdirs.setNoData(FLOWDIR_NO_DATA);

  std::cerr<<"p Calculating D8 flow directions..."<<std::endl;
  progress.start( elevations.width()*elevations.height() );
  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++){
    progress.update( y*elevations.width() );
    for(int x=0;x<elevations.width();x++)
      if(elevations(x,y)==elevations.noData())
        flowdirs(x,y) = flowdirs.noData();
      else
        flowdirs(x,y) = d8_FlowDir(elevations,x,y);
  }
  std::cerr<<"t Succeeded in = "<<progress.stop()<<" s"<<std::endl;
}

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";
}

void dinf_flow_directions(const Array2D<T> &elevations, Array2D<float> &flowdirs){
  ProgressBar progress;

  std::cerr<<"\nA Dinf Flow Directions"<<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 Dinf flow directions matrix..."<<std::endl;
  flowdirs.resize(elevations);
  flowdirs.setNoData(dinf_NO_DATA);
  flowdirs.setAll(NO_FLOW);

  std::cerr<<"p Calculating Dinf flow directions..."<<std::endl;
  progress.start( elevations.size() );
  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++){
    progress.update( y*elevations.width() );
    for(int x=0;x<elevations.width();x++)
      if(elevations(x,y)==elevations.noData())
        flowdirs(x,y) = flowdirs.noData();
      else
        flowdirs(x,y) = dinf_FlowDir(elevations,x,y);
  }
  std::cerr<<"t Succeeded in = "<<progress.stop()<<" s"<<std::endl;
}

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 priority_flood_watersheds(
    Array2D<elev_t> &elevations, Array2D<int32_t> &labels, bool alter_elevations
)
{
    grid_cellz_pq<elev_t> open;
    std::queue<grid_cellz<elev_t> > pit;
    unsigned long processed_cells = 0;
    unsigned long pitc = 0, openc = 0;
    int clabel = 1;  //TODO: Thought this was more clear than zero in the results.
    ProgressBar progress;

    std::cerr << "\n###Priority-Flood+Watershed Labels" << std::endl;
    std::cerr << "Setting up boolean flood array matrix..." << std::flush;
    Array2D<int8_t> closed(elevations.viewWidth(), elevations.viewHeight(), false);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "Setting up watershed label matrix..." << std::flush;
    labels.resize(elevations.viewWidth(), elevations.viewHeight(), -1);
    labels.setNoData(-1);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "The priority queue will require approximately "
              << (elevations.viewWidth() * 2 + elevations.viewHeight() * 2)*((long)sizeof(grid_cellz<elev_t>)) / 1024 / 1024
              << "MB of RAM."
              << std::endl;
    std::cerr << "Adding cells to the priority queue..." << std::endl;
    for (int x = 0; x < elevations.viewWidth(); x++)
    {
        open.push_cell(x, 0, elevations(x, 0));
        open.push_cell(x, elevations.viewHeight() - 1, elevations(x, elevations.viewHeight() - 1));
        closed(x, 0) = true;
        closed(x, elevations.viewHeight() - 1) = true;
    }
    for (int y = 1; y < elevations.viewHeight() - 1; y++)
    {
        open.push_cell(0, y, elevations(0, y));
        open.push_cell(elevations.viewWidth() - 1, y, elevations(elevations.viewWidth() - 1, y));
        closed(0, y) = true;
        closed(elevations.viewWidth() - 1, y) = true;
    }
    std::cerr << "succeeded." << std::endl;

    std::cerr << "%%Performing Priority-Flood+Watershed Labels..." << std::endl;
    progress.start(elevations.viewWidth()*elevations.viewHeight());
    while (open.size() > 0 || pit.size()>0)
    {
        grid_cellz<elev_t> c;
        if (pit.size() > 0)
        {
            c = pit.front();
            pit.pop();
            pitc++;
        }
        else
        {
            c = open.top();
            open.pop();
            openc++;
        }
        processed_cells++;

        //Since all interior cells will be flowing into a cell which has already
        //been processed, the following line identifies only the edge cells of the
        //DEM. Each edge cell seeds its own watershed/basin. The result of this will
        //be many small watersheds/basins around the edge of the DEM.
        if (labels(c.x, c.y) == labels.noData() && elevations(c.x, c.y) != elevations.noData())  //Implies a cell without a label which borders the edge of the DEM or a region of no_data
            labels(c.x, c.y) = clabel++;

        for (int n = 1; n <= 8; n++)
        {
            int nx = c.x + dx[n];
            int ny = c.y + dy[n];
            if (!elevations.in_grid(nx, ny)) continue;
            if (closed(nx, ny))
                continue;

            //Since the neighbouring cell is not closed, its flow is directed to this
            //cell. Therefore, it is part of the same watershed/basin as this cell.
            labels(nx, ny) = labels(c.x, c.y);

            closed(nx, ny) = true;
            if (elevations(nx, ny) <= c.z)
            {
                if (alter_elevations)
                    elevations(nx, ny) = c.z;
                pit.push(grid_cellz<elev_t>(nx, ny, c.z));
            }
            else
                open.push(grid_cellz<elev_t>(nx, ny, elevations(nx, ny)));
        }
        progress.update(processed_cells);
    }

    std::cerr << "\t\033[96msucceeded in " << progress.stop() << "s.\033[39m" << std::endl;

    std::cerr << processed_cells << " cells processed. "
              << pitc << " in pits, "
              << openc << " not in pits."
              << std::endl;
}

void pit_mask(const Array2D<elev_t> &elevations, Array2D<int32_t> &pit_mask)
{
    grid_cellz_pq<elev_t> open;
    std::queue<grid_cellz<elev_t> > pit;
    unsigned long processed_cells = 0;
    unsigned long pitc = 0;
    ProgressBar progress;

    std::cerr << "\n###Pit Mask" << std::endl;
    std::cerr << "Setting up boolean flood array matrix..." << std::flush;
    Array2D<int8_t> closed(elevations.viewWidth(), elevations.viewHeight(), false);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "Setting up the pit mask matrix..." << std::endl;
    pit_mask.resize(elevations.viewWidth(), elevations.viewHeight());
    pit_mask.setNoData(3);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "The priority queue will require approximately "
              << (elevations.viewWidth() * 2 + elevations.viewHeight() * 2)*((long)sizeof(grid_cellz<elev_t>)) / 1024 / 1024
              << "MB of RAM."
              << std::endl;
    std::cerr << "Adding cells to the priority queue..." << std::flush;
    for (int x = 0; x < elevations.viewWidth(); x++)
    {
        open.push_cell(x, 0, elevations(x, 0));
        open.push_cell(x, elevations.viewHeight() - 1, elevations(x, elevations.viewHeight() - 1));
        closed(x, 0) = true;
        closed(x, elevations.viewHeight() - 1) = true;
    }
    for (int y = 1; y < elevations.viewHeight() - 1; y++)
    {
        open.push_cell(0, y, elevations(0, y));
        open.push_cell(elevations.viewWidth() - 1, y, elevations(elevations.viewWidth() - 1, y));
        closed(0, y) = true;
        closed(elevations.viewWidth() - 1, y) = true;
    }
    std::cerr << "succeeded." << std::endl;

    std::cerr << "%%Performing the pit mask..." << std::endl;
    progress.start(elevations.viewWidth()*elevations.viewHeight());
    while (open.size() > 0 || pit.size()>0)
    {
        grid_cellz<elev_t> c;
        if (pit.size() > 0)
        {
            c = pit.front();
            pit.pop();
        }
        else
        {
            c = open.top();
            open.pop();
        }
        processed_cells++;

        for (int n = 1; n <= 8; n++)
        {
            int nx = c.x + dx[n];
            int ny = c.y + dy[n];
            if (!elevations.in_grid(nx, ny)) continue;
            if (closed(nx, ny))
                continue;

            closed(nx, ny) = true;
            if (elevations(nx, ny) <= c.z)
            {
                if (elevations(nx, ny) < c.z)
                    pit_mask(nx, ny) = 1;
                pit.push(grid_cellz<elev_t>(nx, ny, c.z));
            }
            else
            {
                pit_mask(nx, ny) = 0;
                open.push_cell(nx, ny, elevations(nx, ny));
            }
        }

        if (elevations(c.x, c.y) == elevations.noData())
            pit_mask(c.x, c.y) = pit_mask.noData();

        progress.update(processed_cells);
    }
    std::cerr << "\t\033[96msucceeded in " << progress.stop() << "s.\033[39m" << std::endl;
    std::cerr << processed_cells << " cells processed. " << pitc << " in pits." << std::endl;
}

void priority_flood_flowdirs(const Array2D<elev_t> &elevations, Array2D<int8_t> &flowdirs)
{
    grid_cellzk_pq<elev_t> open;
    unsigned long processed_cells = 0;
    ProgressBar progress;

    std::cerr << "\n###Priority-Flood+Flow Directions" << std::endl;
    std::cerr << "Setting up boolean flood array matrix..." << std::flush;
    Array2D<int8_t> closed(elevations.viewWidth(), elevations.viewHeight(), false);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "Setting up the flowdirs matrix..." << std::flush;
    flowdirs.resize(elevations.viewWidth(), elevations.viewHeight());
    flowdirs.setNoData(NO_FLOW);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "The priority queue will require approximately "
              << (elevations.viewWidth() * 2 + elevations.viewHeight() * 2)*((long)sizeof(grid_cellz<elev_t>)) / 1024 / 1024
              << "MB of RAM."
              << std::endl;

    std::cerr << "Adding cells to the priority queue..." << std::endl;
    for (int x = 0; x < elevations.viewWidth(); x++)
    {
        open.push_cell(x, 0, elevations(x, 0));
        open.push_cell(x, elevations.viewHeight() - 1, elevations(x, elevations.viewHeight() - 1));
        flowdirs(x, 0) = 3;
        flowdirs(x, elevations.viewHeight() - 1) = 7;
        closed(x, 0) = true;
        closed(x, elevations.viewHeight() - 1) = true;
    }
    for (int y = 1; y < elevations.viewHeight() - 1; y++)
    {
        open.push_cell(0, y, elevations(0, y));
        open.push_cell(elevations.viewWidth() - 1, y, elevations(elevations.viewWidth() - 1, y));
        flowdirs(0, y) = 1;
        flowdirs(elevations.viewWidth() - 1, y) = 5;
        closed(0, y) = true;
        closed(elevations.viewWidth() - 1, y) = true;
    }
    std::cerr << "succeeded." << std::endl;

    flowdirs(0, 0) = 2;
    flowdirs(flowdirs.viewWidth() - 1, 0) = 4;
    flowdirs(0, flowdirs.viewHeight() - 1) = 8;
    flowdirs(flowdirs.viewWidth() - 1, flowdirs.viewHeight() - 1) = 6;

    const int d8_order[9] = { 0,1,3,5,7,2,4,6,8 };
    std::cerr << "%%Performing Priority-Flood+Flow Directions..." << std::endl;
    progress.start(elevations.viewWidth()*elevations.viewHeight());
    while (open.size() > 0)
    {
        grid_cellz<elev_t> c = open.top();
        open.pop();
        processed_cells++;

        for (int no = 1; no <= 8; no++)
        {
            int n = d8_order[no];
            int nx = c.x + dx[n];
            int ny = c.y + dy[n];
            if (!elevations.in_grid(nx, ny)) continue;
            if (closed(nx, ny))
                continue;

            closed(nx, ny) = true;

            if (elevations(nx, ny) == elevations.noData())
                flowdirs(nx, ny) = flowdirs.noData();
            else
                flowdirs(nx, ny) = inverse_flow[n];

            open.push_cell(nx, ny, elevations(nx, ny));
        }
        progress.update(processed_cells);
    }
    std::cerr << "\t\033[96msucceeded in " << progress.stop() << "s.\033[39m" << std::endl;
    std::cerr << processed_cells << " cells processed." << std::endl;
}