What is a simple binary raster storage format which ArcGIS can read?

2013-05-10

An answer to this question on the GIS Stack Exchange.

Question

I have a C++ program that manipulates raster digital elevation models.

I have been using ArcGrid ASCII format for results. It is simple to code and gets the job done, but it necessitates a decimal output. This means a loss of precision in floating-point numbers.

Is there a simple binary raster format which I could write to and ArcGIS could read from?

Answer

I ended up using the ArcGIS floating-point raster format (extension .flt), as described here).

The C++ code I wrote to handle things follows. It will not compile without a grid class, but it should be easy enough to read and understand.

/**
  @brief  Writes a floating-point grid file
  @param[in]  &basename     Name, without extension, of output file
  @param[in]  &output_grid  DEM object to write
  @todo Does not check byte order (big-endian, little-endian)
  @todo Does not output only IEEE-754 32-bit floating-point,
        which is required by ArcGIS
  @returns 0 upon success
*/
template <class T>
int write_floating_data(const std::string basename,const array2d<T> &output_grid){
  std::string fn_header(basename), fn_data(basename);
  fn_header+=".hdr";
  fn_data+=".flt";
  {
    std::ofstream fout;
    fout.open(fn_header.c_str());
    if(!fout.is_open())
      exit(-1);
    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: Detect endian-ness
    fout.close();
  }
  {
    std::ofstream fout(fn_data.c_str(), std::ios::binary | std::ios::out);
    if(!fout.is_open())
      exit(-1);
    for(int y=0;y<output_grid.height();++y)
    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();
  }
  return 0;
}
/**
  @brief  Reads in a floating-point grid file
  @param[in]  &basename     Name, without extension, of input file
  @param[in]  &grid         DEM object in which to store data
  @todo Does not check byte order (big-endian, little-endian)
  @todo Does not input only IEEE-754 32-bit floating-point,
        which is required by ArcGIS
  @returns 0 upon success
*/
template <class T>
int read_floating_data(const std::string basename,array2d<T> &grid){
  std::string fn_header(basename), fn_data(basename);
  fn_header+=".hdr";
  fn_data+=".flt";
  int columns, rows;
  char byteorder; //L=LSB first, M=MSB first
  {
    FILE *fin;
    fin=fopen(fn_header.c_str(),"r");
    if(fin==NULL)
      exit(-1);
    if(fscanf(fin,"ncols %d nrows %d xllcorner %lf yllcorner %lf cellsize %lf NODATA_value %f BYTEORDER %c",&columns, &rows, &grid.xllcorner, &grid.yllcorner, &grid.cellsize, &grid.no_data, &byteorder)!=7)
      exit(-1);
    fclose(fin);
  }
  grid.resize(columns,rows);
  {
    std::ifstream fin(fn_data.c_str(), std::ios::binary | std::ios::in);
    if(!fin.is_open())
      exit(-1);
    grid.data_cells=0;
    for(int y=0;y<rows;++y)
    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++;
    }
  }
  return 0;
}