C++｜基于D8和LCP algorithm的径流计算实现

Runoff modelling

GitHub

git@github.com:YitongXia/runoff_modelling.git

环境配置

需要安装GDAL
（VS2019可使用vcpkg进行安装）

算法思路

利用D8径流模型（single flow direction, SFD）和LCP algorithm进行计算flow direction，根据生成的direction进行流量汇总，生成flow accumulation

code implementation

include

#include <iostream>
#include <queue>
#include <stack>
#include <cassert>
#include <fstream>
#include "gdal_priv.h"
#include "cpl_conv.h"

struct

( Raster and RasterCell.
Raster is using for storing all the pixel values in the tiff file, RasterCell is for storing information of each pixel)

// A structure that links to a single cell in a Raster
struct RasterCell {
   
    int x, y; // row and column of the cell
    int elevation;
    int insertion_order;
    bool add_to_list;
    int flow_dir;
    int accumulation = 0;

    // Defines a new link to a cell
    RasterCell(int x, int y, int elevation, int insert_order, int dir) {
   
        this->x = x;
        this->y = y;
        this->elevation = elevation;
        this->insertion_order = insert_order;
        this->flow_dir = dir;
    }

    // Define the order of the linked cells (to be used in a priority_queue)
    bool operator<(const RasterCell& other) const {
   
        // to do with statements like if (this->elevation > other.elevation) return false/true;
        if (this->elevation == other.elevation)
        {
   
            if (this->insertion_order < other.insertion_order)
                return false;
            else if (this->insertion_order > other.insertion_order)
                return true;
        }
        else if (this->elevation <= other.elevation)
            return false;
        else if (this->elevation > other.elevation)
            return true;
       
    }
};

// Storage and access of a raster of a given size
struct Raster {
   
    std::vector<unsigned int> pixels; // where everything is stored
    std::vector<int> visiting; // store information about if the cell is visited
    std::vector<int> in_queue; // store information about if the cell is added to the priority queue


    int max_x, max_y; // number of columns and rows
    int direction;
    

    // Initialise a raster with x columns and y rows
    Raster(int x, int y) {
   
        max_x = x;
        max_y = y;
        unsigned int total_pixels = x * y;
        pixels.reserve(total_pixels);
        visiting.reserve(total_pixels);
        
    }

    // Fill values of an entire row
    void add_scanline(const unsigned int* line) {
   
        for (int i = 0; i < max_x; ++i)
            pixels.push_back(line[i]);
    }

    // Fill entire raster with zeros
    void fill() {
   
        unsigned int total_pixels = max_x * max_y;
        for (int i = 0; i < total_pixels; ++i)
        {
   
            pixels.push_back(0);
        }
    }

    // Fill the entire raster with 0 in visiting and in_queue vector
    void fill_visit() {
   
        unsigned int total_pixels = max_x * max_y;
        for (int i = 0; i < total_pixels; ++i)
        {
   
            visiting.push_back(0);
            in_queue.push_back(0);
        }
    }


    // Access the value of a raster cell to read or write it
    unsigned int& operator()(int x, int y) {
   
        assert(x >= 0 && x < max_x);
        assert(y >= 0 && y < max_y);
        return pixels[x + y * max_x];
    }

    // Access the value of a raster cell to read it
    unsigned int operator()(int x, int y) const {
   
        assert(x >= 0 && x < max_x);
        assert(y >= 0 && y < max_y);
        return pixels[x + y * max_x];
    }

    // Add pixel value for output raster
    void add_value(int x1, int y1, int value) {
   
        assert(x1 >= 0 && x1 < max_x);
        assert(y1 >= 0 && y1 < max_y);
        pixels[x1 + y1 * max_x] = value;
        
    }

    //change the status from unvisited to visited
    void Is_Visited(int x1, int y1) {
   
        visiting[x1 + y1 * max_x] = 1;
    }

    // return the status of the pixel whether visited.
    int If_Visited(int x1, int y1) {
   
        return visiting[x1 + y1 * max_x];
    }

    // update the status after add to the cell_to_process_flow queue
    void Add_to_queue(int x, int y) {
   
        in_queue[x + y * max_x] = 1;
    }

    // check the status whether add to the cell_to_process_flow queue
    int If_add_to_queue(int x, int y) {
   
        return in_queue[x + y * max_x];
    }

    void output_accumulation(int& current_line, unsigned int* line)
    {
   
        for (int i = 0; i < max_y; ++i) 
            line[i] = pixels[i + current_line * max_y];
    }
};

other functions and global variants

int insert_order = 0;

// Write the values in a linked raster cell (useful for debugging)
std::ostream& operator<<(std::ostream& os, const RasterCell& c) {
   
    os << "{h=" << c.elevation << ", o=" << c.insertion_order << ", x=" << c.x << ", y=" << c.y << "}";
    return os;
};

// write raster result into the tiff file
void output_tiff(std::string filename, Raster input_raster,int max_x, int max_y) {
   

    int nXSize = max_x;
    int nYSize = max_y;
    double geo_transform[6];
    std::string tiffname = filename;
    GDALAllRegister();
    CPLPushErrorHandler(CPLQuietErrorHandler);
    GDALDataset* gribDataset, * geotiffDataset;
    GDALDriver* driverGeotiff;
    GDALRasterBand* geotiffBand;
    
    gribDataset = (GDALDataset*)GDALOpen("N36E076.hgt", GA_ReadOnly);
    gribDataset->GetGeoTransform(geo_transform);

    driverGeotiff = GetGDALDriverManager()->GetDriverByName("GTiff");
    geotiffDataset = driverGeotiff->Create(tiffname.c_str(), nXSize, nYSize, 1, GDT_Int32, NULL);
    geotiffDataset->SetGeoTransform(geo_transform);
    geotiffDataset->SetProjection(gribDataset->GetProjectionRef());
    unsigned int* rowBuff = (unsigned int*)CPLMalloc(sizeof(unsigned int) * nYSize);
    for (int j = 0; j < nYSize; j++) {
   
        for (int i = 0; i < nXSize; i++) {
   
            rowBuff[i] = (unsigned int)input_raster(i, j);
        }
        geotiffDataset->GetRasterBand(1)->RasterIO(GF_Write, 0, j, nYSize, 1, rowBuff, nYSize, 1, GDT_Int32, 0, 0);
    }
    GDALClose(gribDataset);
    GDALClose(geotiffDataset);
    GDALDestroyDriverManager();
}

calculate flow direction

int main(int argc, const char* argv[]) {
   
    // Open dataset
    GDALDataset* input_dataset;
    GDALAllRegister();
    input_dataset = (GDALDataset*)GDALOpen("N36E076.hgt", GA_ReadOnly);
    if (input_dataset == NULL) {
   
        std::cerr << "Couldn't open file" << std::endl;
        return 1;
    }

    // Print dataset info
    double geo_transform[6];
    std::cout << "Driver: " << input_dataset->GetDriver()->GetDescription() << "/" << input_dataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME) << std::endl;;
    std::cout << "Size is " << input_dataset->GetRasterXSize() << "x" << input_dataset->GetRasterYSize() << "x" << input_dataset->GetRasterCount() << std::endl;
    if (input_dataset->GetProjectionRef() != NULL) std::cout << "Projection is '" << input_dataset->GetProjectionRef() << "'" << std::endl;
    if (input_dataset->GetGeoTransform(geo_transform) == CE_None) {
   
        std::cout << "Origin = (" << geo_transform[0] << ", " << geo_transform[3] << ")" << std::endl;
        std::cout << "Pixel Size = (" << geo_transform[1] << ", " << geo_transform[5] << ")" << std::endl;
    }

    // Print Band 1 info
    GDALRasterBand* input_band;
    int nBlockXSize, nBlockYSize;
    int bGotMin, bGotMax;
    double adfMinMax[2];
    input_band = input_dataset->GetRasterBand(1);
    input_band->GetBlockSize(&nBlockXSize, &nBlockYSize);
    std::cout << "Band 1 Block=" << nBlockXSize << "x" << nBlockYSize << " Type=" << GDALGetDataTypeName(input_band->GetRasterDataType()) << " ColorInterp=" << GDALGetColorInterpretationName(input_band->GetColorInterpretation()) << std::endl;
    adfMinMax[0] = input_band->GetMinimum(&bGotMin);
    adfMinMax[1] = input_band->GetMaximum(&bGotMax);
    if (!(bGotMin && bGotMax))
        GDALComputeRasterMinMax((GDALRasterBandH)input_band, TRUE, adfMinMax);
    std::cout << "Min=" << adfMinMax[0] << " Max=" << adfMinMax[1] << std::endl;

    // Read Band 1 line by line
    int nXSize = input_band->GetXSize();
    int nYSize = input_band->GetYSize();
    Raster input_raster(nXSize, nYSize);
    for (int current_scanline