C++ array2d::width方法代码示例

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

  diagnostic("Setting up the flow directions matrix...");
  flowdirs.copyprops(elevations);
  flowdirs.init(NO_FLOW);
  flowdirs.no_data=d8_NO_DATA;
  diagnostic("succeeded.\n");

  diagnostic("%%Calculating D8 flow directions...\n");
  progress.start( elevations.width()*elevations.height() );
  #pragma omp parallel for
  for(int x=0;x<elevations.width();x++){
    progress.update( x*elevations.height() );
    for(int y=0;y<elevations.height();y++)
      if(elevations(x,y)==elevations.no_data)
        flowdirs(x,y)=flowdirs.no_data;
      else
        flowdirs(x,y)=d8_FlowDir(elevations,x,y);
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());
}

    reference operator=(const array2d<float> arr){
      if(arr.width()!=width() || arr.height()!=height())
        throw "Error! Grids were not of the same size!";

      #pragma omp parallel for collapse(2)
      for(int x=0;x<width();x++)
      for(int y=0;y<height();y++)
        operator()(x,y)=arr(x,y);
    }

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.edge_grid(x,y)){
    if(x==0 && y==0)
      return 2;
    else if(x==0 && y==elevations.height()-1)
      return 8;
    else if(x==elevations.width()-1 && y==0)
      return 4;
    else if(x==elevations.width()-1 && y==elevations.height()-1)
      return 6;
    else if(x==0)
      return 1;
    else if(x==elevations.width()-1)
      return 5;
    else if(y==0)
      return 3;
    else if(y==elevations.height()-1)
      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)
    ){
      minimum_elevation=elevations(x+dx[n],y+dy[n]);
      flowdir=n;
    }

  return flowdir;
}

void d8_upslope_area(const array2d<T> &flowdirs, array2d<U> &area){
  char_2d dependency;
  std::queue<grid_cell> sources;
  ProgressBar progress;

  diagnostic("\n###D8 Upslope Area\n");

  diagnostic_arg(
    "The sources queue will require at most approximately %ldMB of RAM.\n",
    flowdirs.width()*flowdirs.height()*((long)sizeof(grid_cell))/1024/1024
  );

  diagnostic("Resizing dependency matrix...");
  dependency.copyprops(flowdirs);
  diagnostic("succeeded.\n");

  diagnostic("Setting up the area matrix...");
  area.copyprops(flowdirs);
  area.init(0);
  area.no_data=d8_NO_DATA;
  diagnostic("succeeded.\n");

  diagnostic("%%Calculating dependency matrix & setting no_data cells...\n");
  progress.start( flowdirs.width()*flowdirs.height() );
  #pragma omp parallel for
  for(int x=0;x<flowdirs.width();x++){
    progress.update( x*flowdirs.height() );
    for(int y=0;y<flowdirs.height();y++){
      dependency(x,y)=0;
      if(flowdirs(x,y)==flowdirs.no_data){
        area(x,y)=area.no_data;
        continue;
      }
      for(int n=1;n<=8;n++)
        if(!flowdirs.in_grid(x+dx[n],y+dy[n]))
          continue;
        else if(flowdirs(x+dx[n],y+dy[n])==NO_FLOW)
          continue;
        else if(flowdirs(x+dx[n],y+dy[n])==flowdirs.no_data)
          continue;
        else if(n==inverse_flow[(int)flowdirs(x+dx[n],y+dy[n])])
          ++dependency(x,y);
    }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  diagnostic("%%Locating source cells...\n");
  progress.start( flowdirs.width()*flowdirs.height() );
  for(int x=0;x<flowdirs.width();x++){
    progress.update( x*flowdirs.height() );
    for(int y=0;y<flowdirs.height();y++)
      if(flowdirs(x,y)==flowdirs.no_data)
        continue;
      else if(dependency(x,y)==0)
        sources.push(grid_cell(x,y));
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  diagnostic("%%Calculating up-slope areas...\n");
  progress.start(flowdirs.data_cells);
  long int ccount=0;
  while(sources.size()>0){
    grid_cell c=sources.front();
    sources.pop();

    ccount++;
    progress.update(ccount);

    area(c.x,c.y)+=1;

    if(flowdirs(c.x,c.y)==NO_FLOW)
      continue;

    int nx=c.x+dx[(int)flowdirs(c.x,c.y)];
    int ny=c.y+dy[(int)flowdirs(c.x,c.y)];
    if(flowdirs.in_grid(nx,ny) && area(nx,ny)!=area.no_data){
      area(nx,ny)+=area(c.x,c.y);
      if((--dependency(nx,ny))==0)
        sources.push(grid_cell(nx,ny));
    }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());
}

int load_ascii_data(std::string filename, array2d<T> &elevations){
  std::ifstream fin;
  size_t file_size;
  int rows,columns;
  Timer load_time;
  ProgressBar progress;

  load_time.start();

  diagnostic_arg("Opening input ASCII-DEM file \"%s\"...",filename.c_str());
  fin.open(filename.c_str());
  if(!fin.good()){
    diagnostic("failed!\n");
    exit(-1);
  }
  diagnostic("succeeded.\n");

  diagnostic("Calculating file size...");
  fin.seekg(0, fin.end);
  file_size=fin.tellg();
  fin.seekg(0, fin.beg);
  diagnostic("succeeded.\n");

//  posix_fadvise(fileno(fin),0,0,POSIX_FADV_SEQUENTIAL);

  diagnostic("Reading DEM header...");
  fin>>must_be("ncols")         >>columns;
  fin>>must_be("nrows")         >>rows;
  fin>>must_be("xllcorner")     >>elevations.xllcorner;
  fin>>must_be("yllcorner")     >>elevations.yllcorner;
  fin>>must_be("cellsize")      >>elevations.cellsize;
  try {
    fin>>must_be("NODATA_value")  >>elevations.no_data;
  } catch (std::string e) {
    std::cerr<<e<<std::endl;
    std::cerr<<"Continuing without a NoData value!"<<std::endl;
  }
  diagnostic("succeeded.\n");

  diagnostic_arg("The loaded DEM will require approximately %ldMB of RAM.\n",columns*rows*((long)sizeof(float))/1024/1024);

  diagnostic("Resizing elevation matrix...");  //TODO: Consider abstracting this block
  elevations.resize(columns,rows);
  diagnostic("succeeded.\n");

  diagnostic("%%Reading elevation matrix...\n");
  progress.start(file_size);

  elevations.data_cells=0;
  for(int y=0;y<rows;y++){
    progress.update(fin.tellg()); //Todo: Check to see if ftell fails here?
    for(int x=0;x<columns;x++){
      fin>>elevations(x,y);
      if(elevations(x,y)!=elevations.no_data)
        elevations.data_cells++;
    }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  fin.close();

  diagnostic_arg(
    "Read %ld cells, of which %ld contained data (%ld%%).\n",
    elevations.width()*elevations.height(), elevations.data_cells,
    elevations.data_cells*100/elevations.width()/elevations.height()
  );

  load_time.stop();
  diagnostic_arg("Read time was: %lfs\n", load_time.accumulated());

  return 0;
}

int read_floating_data(
  const std::string basename,
  array2d<T> &grid
){
  Timer io_time;
  ProgressBar progress;
  std::string fn_header(basename), fn_data(basename);

  //TODO: The section below should work, but is something of an abomination
  if(typeid(T)==typeid(float)){
    fn_header+=".hdr";
    fn_data+=".flt";
  } else if (typeid(T)==typeid(double)){
    fn_header+=".hdr";
    fn_data+=".dflt";
  } else {
    std::cerr<<"Cannot read floating type data into this format!"<<std::endl;
    exit(-1);
  }

  int columns, rows;
  std::string byteorder;

  io_time.start();


  {
    std::ifstream fin;
    diagnostic_arg("Opening floating-point header file \"%s\" for reading...",fn_header.c_str());
    fin.open(fn_header.c_str());
    if(fin==NULL){
      diagnostic("failed!\n");
      exit(-1);
    }
    diagnostic("succeeded.\n");


    diagnostic("Reading DEM header...");
    fin>>must_be("ncols")         >>columns;
    fin>>must_be("nrows")         >>rows;
    fin>>must_be("xllcorner")     >>grid.xllcorner;
    fin>>must_be("yllcorner")     >>grid.yllcorner;
    fin>>must_be("cellsize")      >>grid.cellsize;
    fin>>must_be("NODATA_value")  >>grid.no_data;
    fin>>must_be("BYTEORDER")     >>byteorder;
    diagnostic("succeeded.\n");
    fin.close();
  }

  diagnostic_arg("The loaded DEM will require approximately %ldMB of RAM.\n",columns*rows*((long)sizeof(float))/1024/1024);

  diagnostic("Resizing grid...");  //TODO: Consider abstracting this block
  grid.resize(columns,rows);
  diagnostic("succeeded.\n");



  diagnostic_arg("Opening floating-point data file \"%s\" for reading...",fn_data.c_str());

  {
    std::ifstream fin(fn_data.c_str(), std::ios::binary | std::ios::in);
    if(!fin.is_open()){
      diagnostic("failed!\n");
      exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
    }
    diagnostic("succeeded.\n");


    diagnostic("%%Reading data...\n");
    progress.start(columns*rows);
    grid.data_cells=0;
    for(int y=0;y<rows;++y){
      progress.update(y*columns); //Todo: Check to see if ftell fails here?
      for(int x=0;x<columns;++x){
        fin.read(reinterpret_cast<char*>(&grid(x,y)), std::streamsize(sizeof(T)));
        if(grid(x,y)!=grid.no_data)
          grid.data_cells++;
      }
    }
    io_time.stop();
    diagnostic_arg(SUCCEEDED_IN,progress.stop());

  }

  diagnostic_arg(
    "Read %ld cells, of which %ld contained data (%ld%%).\n",
    grid.width()*grid.height(), grid.data_cells,
    grid.data_cells*100/grid.width()/grid.height()
  );

  diagnostic_arg("Read time was: %lf\n", io_time.accumulated());

  return 0;
}

int write_floating_data(
  const std::string basename,
  const array2d<T> &output_grid
){
  Timer write_time;
  ProgressBar progress;
  std::string fn_header(basename), fn_data(basename);

  //TODO: The section below should work, but is something of an abomination
  if(typeid(T)==typeid(float)){
    fn_header+=".hdr";
    fn_data+=".flt";
  } else if (typeid(T)==typeid(double)){
    fn_header+=".hdr";
    fn_data+=".dflt";
  } else {
    std::cerr<<"Cannot write floating type data into this format!"<<std::endl;
    exit(-1);
  }

  write_time.start();


  {
    diagnostic_arg("Opening floating-point header file \"%s\" for writing...",fn_header.c_str());
    std::ofstream fout;
    fout.open(fn_header.c_str());
    if(!fout.is_open()){
      diagnostic("failed!\n");
      exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
    }
    diagnostic("succeeded.\n");

    diagnostic("Writing floating-point header file...");
    fout<<"ncols\t\t"<<output_grid.width()<<std::endl;
    fout<<"nrows\t\t"<<output_grid.height()<<std::endl;
    fout<<"xllcorner\t"<<std::fixed<<std::setprecision(10)<<output_grid.xllcorner<<std::endl;
    fout<<"yllcorner\t"<<std::fixed<<std::setprecision(10)<<output_grid.yllcorner<<std::endl;
    fout<<"cellsize\t"<<std::fixed<<std::setprecision(10)<<output_grid.cellsize<<std::endl;
    fout<<"NODATA_value\t"<<std::fixed<<std::setprecision(10)<<output_grid.no_data<<std::endl;
    fout<<"BYTEORDER\tLSBFIRST"<<std::endl; //TODO
    fout.close();
    diagnostic("succeeded.\n");
  }


  diagnostic_arg("Opening floating-point data file \"%s\" for writing...",fn_data.c_str());

  {
    std::ofstream fout(fn_data.c_str(), std::ios::binary | std::ios::out);
    if(!fout.is_open()){
      diagnostic("failed!\n");
      exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
    }
    diagnostic("succeeded.\n");

    diagnostic("%%Writing floating-point data file...\n");
    progress.start( output_grid.width()*output_grid.height() );
    for(int y=0;y<output_grid.height();++y){
      progress.update( y*output_grid.width() );
      for(int x=0;x<output_grid.width();++x)
        fout.write(reinterpret_cast<const char*>(&output_grid(x,y)), std::streamsize(sizeof(T)));
    }
    fout.close();
    write_time.stop();
    diagnostic_arg(SUCCEEDED_IN,progress.stop());
  }

  diagnostic_arg("Write time was: %lf\n", write_time.accumulated());

  return 0;
}

int output_ascii_data(
  const std::string filename,
  const array2d<T> &output_grid,
  int precision=8
){
  std::ofstream fout;
  std::string outputsep=" ";
  int output_type=OUTPUT_DEM;
  Timer write_time;
  ProgressBar progress;

  write_time.start();

  diagnostic_arg("Opening ASCII output file \"%s\"...",filename.c_str());
  fout.open(filename.c_str());
  if(!fout.is_open()){
    diagnostic("failed!\n");
    exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
  }
  diagnostic("succeeded.\n");

  //OmniGlyph output
  if(filename.substr(filename.length()-4)==".omg"){
    outputsep="|";
    output_type=OUTPUT_OMG;
    diagnostic("Writing OmniGlyph file header...");
    fout<<"Contents: Pixel array"<<std::endl;
    fout<<std::endl;
    fout<<"Width:    "<<output_grid.width()<<std::endl;
    fout<<"Height:   "<<output_grid.height()<<std::endl;
    fout<<std::endl;
    fout<<"Spectral bands:   1"<<std::endl;
    fout<<"Bits per band:   32"<<std::endl;
    fout<<"Range of values:   "<<output_grid.min()<<","<<output_grid.max()<<std::endl;
    fout<<"Actual range:   "<<output_grid.no_data<<","<<output_grid.max()<<std::endl;  //TODO: Assumes no_data is a small negative value
    fout<<"Gamma exponent:   0."<<std::endl;
    fout<<"Resolution:   100 pixels per inch"<<std::endl;
    fout<<std::endl;
    fout<<"|"<<std::endl;
  } else {
    diagnostic("Writing ArcGrid ASCII file header...");
    fout<<"ncols\t\t"<<output_grid.width()<<std::endl;
    fout<<"nrows\t\t"<<output_grid.height()<<std::endl;
    fout<<"xllcorner\t"<<std::fixed<<std::setprecision(precision)<<output_grid.xllcorner<<std::endl;
    fout<<"yllcorner\t"<<std::fixed<<std::setprecision(precision)<<output_grid.yllcorner<<std::endl;
    fout<<"cellsize\t"<<std::fixed<<std::setprecision(precision)<<output_grid.cellsize<<std::endl;
    fout<<"NODATA_value\t"<<std::fixed<<std::setprecision(precision);
    if(sizeof(T)==1)  //TODO: Crude way of detecting chars and bools
      fout<<(int)output_grid.no_data<<std::endl;
    else
      fout<<output_grid.no_data<<std::endl;
  }
  diagnostic("succeeded.\n");

  diagnostic("%%Writing ArcGrid ASCII file data...\n");
  progress.start( output_grid.width()*output_grid.height() );
  fout.precision(precision);
  fout.setf(std::ios::fixed);
  for(int y=0;y<output_grid.height();y++){
    progress.update( y*output_grid.width() );
    if(output_type==OUTPUT_OMG)
      fout<<"|";
    for(int x=0;x<output_grid.width();x++)
      if(sizeof(T)==1)  //TODO: Crude way of detecting chars and bools
        fout<<(int)output_grid(x,y)<<outputsep;
      else
        fout<<output_grid(x,y)<<outputsep;
    fout<<std::endl;
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  fout.close();

  write_time.stop();
  diagnostic_arg("Write time was: %lf\n", write_time.accumulated());

  return 0;
}