【QGIS二次开发】空间分析-10_qgis 示例代码-优快云博客

系列目录：

5.1 矢量缓冲区分析

功能点描述：

矢量缓冲区分析：分别实现点、线路和区的缓冲区，并输出结果区文件。

运行截图：

操作前：

图 89 示例数据

操作后：

图 90 缓冲区结果

主要代码如下：

该代码是一个执行缓冲区分析的函数performBufferAnalysis。首先，它通过GDAL库注册所有驱动器，打开输入文件，获取输入图层，并创建输出文件和输出图层。接着，它使用OGRGeometryFactory创建一个用于缓冲区分析的几何对象(wkbPolygon)。然后，通过遍历输入图层的要素，获取每个要素的几何对象，并对几何对象执行缓冲区分析。缓冲距离和段数通过参数传递，生成缓冲区分析后的几何对象。

随后，创建输出要素，将缓冲区分析后的几何对象设置为输出要素的几何对象，并将输出要素写入输出图层。最后，释放内存，关闭输入文件和输出文件，弹出提示框告知用户缓冲区分析已完成。然后，通过QgsVectorLayer类将输出文件作为图层添加到QgsProject中，并刷新地图画布，以在地图中显示缓冲区分析的结果。

void ZoneAnalysis::performBufferAnalysis(const QString& inputFilePath, const QString& outputFilePath, double bufferDistance) {    
    // 注册所有的驱动器    
    GDALAllRegister();    
    
    // 打开输入文件    
    GDALDataset* inputDataset = (GDALDataset*)GDALOpenEx(inputFilePath.toStdString().c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr);    
    if (inputDataset == nullptr) {    
        QMessageBox::critical(this, "错误", "无法打开输入文件！");    
        return;    
    }    
    
    // 获取输入图层    
    OGRLayer* inputLayer = inputDataset->GetLayer(0);    
    
    // 创建输出文件    
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");    
    GDALDataset* outputDataset = driver->Create(outputFilePath.toStdString().c_str(), 0, 0, 0, GDT_Unknown, nullptr);    
    if (outputDataset == nullptr) {    
        QMessageBox::critical(this, "错误", "无法创建输出文件！");    
        GDALClose(inputDataset);    
        return;    
    }    
    
    
    // 创建输出图层    
    OGRLayer* outputLayer = outputDataset->CreateLayer("buffer", nullptr, wkbPolygon, nullptr);    
    if (outputLayer == nullptr) {    
        QMessageBox::critical(this, "错误", "无法创建输出图层！");    
        GDALClose(outputDataset);    
        GDALClose(inputDataset);    
        return;    
    }    
    
    // 创建缓冲区分析对象    
    OGRGeometryFactory::createGeometry(wkbPolygon);    
    
    // 遍历输入图层的要素    
    inputLayer->ResetReading();    
    OGREnvelope envelope;    
    while (OGRFeature* feature = inputLayer->GetNextFeature()) {    
        OGRGeometry* geometry = feature->GetGeometryRef();    
        if (geometry != nullptr) {    
            // 执行缓冲区分析    
            OGRGeometry* bufferGeometry = geometry->Buffer(bufferDistance, 30);    
    
            // 创建输出要素    
            OGRFeature* outputFeature = OGRFeature::CreateFeature(outputLayer->GetLayerDefn());    
            outputFeature->SetGeometry(bufferGeometry);    
    
            // 将输出要素写入输出图层    
            outputLayer->CreateFeature(outputFeature);    
    
            // 释放内存    
            OGRFeature::DestroyFeature(outputFeature);    
            OGRGeometryFactory::destroyGeometry(bufferGeometry);    
        }    
    
        // 释放内存    
        OGRFeature::DestroyFeature(feature);    
    }    
    
    // 释放资源    
    GDALClose(outputDataset);    
    GDALClose(inputDataset);    
    // 缓冲区分析完成后，关闭缓冲区框    
    this->close();    
    
    // 弹出提示框告知用户    
    QMessageBox::information(this, "提示", "缓冲区分析已完成！");    
    
    // 从输出路径打开shp    
    QFileInfo fi(outputFilePath);    
    if (!fi.exists()) { return; }    
    QString layerBaseName = fi.baseName(); // 图层名称    
    QgsVectorLayer* vecLayer = new QgsVectorLayer(outputFilePath, layerBaseName);    
    if (!vecLayer) { return; }    
    if (!vecLayer->isValid())    
    {    
        QMessageBox::information(0, "", "layer is invalid");    
        return;    
    }    
    QgsProject::instance()->addMapLayer(vecLayer);    
    m_mapCanvas->refresh();    
    //m_mapCanvas->setExtent(vecLayer->extent());    
    //layers.append(vecLayer);    
    //m_mapCanvas->setLayers(layers);    
    //m_mapCanvas->refresh();    
}

5.2 矢量裁剪分析

功能点描述：

矢量裁剪分析：分别用矩形、圆和任意多边形对矢量图层进行裁剪，并将裁剪结果输出到矢量文件。

运行截图：

操作前：

图 91 裁剪数据

操作后：

图 92 裁剪效果

主要代码：

这段代码是一个名为ZoneCut的类的成员函数，用于对两个矢量文件进行裁剪，并将裁剪后的结果输出到一个新的矢量文件中。首先，它初始化GDAL库。然后，通过打开两个输入文件，创建输出文件，以及获取相应的图层等操作，准备进行矢量文件的裁剪操作。如果打开输入文件或创建输出文件失败，则函数会提前返回。

在裁剪操作的过程中，它遍历第一个输入文件的每一个特征，获取其几何体，并为每个特征创建一个裁剪几何体的克隆。接着，将第二个输入文件的空间过滤器设置为裁剪几何体，遍历第二个输入文件的每一个特征，如果其几何体与裁剪几何体相交，则计算它们的交集，并将交集作为输出图层的一个特征。

最后，完成裁剪操作后，关闭输入文件和输出文件，清空输入和输出文件路径，同时关闭与裁剪相关的对话框。该函数主要实现了两个矢量文件的裁剪，将裁剪后的结果保存到一个新的矢量文件中，为用户提供了对地理数据进行空间裁剪的功能。

void ZoneCut::performZoneCut(const QString& inputFilePath1, const QString& inputFilePath2, const QString& outputFilePath)    
{    
    GDALAllRegister();    
    // 打开第一个输入文件    
    QString inputFilePath1Utf8 = inputFilePath1.toUtf8();    
    GDALDataset* inputDataset1 = (GDALDataset*)GDALOpenEx(inputFilePath1.toStdString().c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr);    
    if (inputDataset1 == nullptr)    
    {    
        // 处理打开文件失败的情况    
        return;    
    }    
    
    // 打开第二个输入文件    
    GDALDataset* inputDataset2 = (GDALDataset*)GDALOpenEx(inputFilePath2.toStdString().c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr);    
    if (inputDataset2 == nullptr)    
    {    
        // 处理打开文件失败的情况    
        GDALClose(inputDataset1);    
        return;    
    }    
    
    // 创建输出文件    
    GDALDriver* outputDriver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");    
    GDALDataset* outputDataset = outputDriver->Create(outputFilePath.toStdString().c_str(), 0, 0, 0, GDT_Unknown, nullptr);    
    if (outputDataset == nullptr)    
    {    
        // 处理创建文件失败的情况    
        GDALClose(inputDataset1);    
        GDALClose(inputDataset2);    
        return;    
    }    
    
    // 获取输入文件的第一个图层    
    OGRLayer* inputLayer1 = inputDataset1->GetLayer(0);    
    
    // 获取输入文件的第二个图层    
    OGRLayer* inputLayer2 = inputDataset2->GetLayer(0);    
    
    // 创建输出图层    
    OGRLayer* outputLayer = outputDataset->CreateLayer("result", inputLayer1->GetSpatialRef()->Clone(), inputLayer1->GetLayerDefn()->GetGeomType(), nullptr);    
    
    // 获取第一个输入图层的几何体    
    OGRFeature* inputFeature1 = nullptr;//初始化输入文件1的特征指针。    
    while ((inputFeature1 = inputLayer1->GetNextFeature()) != nullptr)    
    {    
        OGRGeometry* inputGeometry1 = inputFeature1->GetGeometryRef();    
        OGRGeometry* cutGeometry = inputGeometry1->clone();//克隆输入文件1的特征的几何体    
    
        // 设置第二个输入图层的空间过滤器为被裁剪几何体    
        inputLayer2->SetSpatialFilter(cutGeometry);    
    
        // 获取第二个输入图层的几何体    
        OGRFeature* inputFeature2 = nullptr;    
        while ((inputFeature2 = inputLayer2->GetNextFeature()) != nullptr)    
        {    
            OGRGeometry* inputGeometry2 = inputFeature2->GetGeometryRef();// 获取输入文件2的特征的几何体。    
            if (inputGeometry2->Intersects(cutGeometry))//如果输入文件2的特征的几何体与裁剪几何体相交。    
            {    
                OGRFeature* outputFeature = OGRFeature::CreateFeature(outputLayer->GetLayerDefn());    
                outputFeature->SetGeometry(inputGeometry2->Intersection(cutGeometry));//设置输出图层的特征的几何体为输入文件2的特征的几何体与裁剪几何体的交集。    
                outputLayer->CreateFeature(outputFeature);    
                OGRFeature::DestroyFeature(outputFeature);    
            }    
            OGRFeature::DestroyFeature(inputFeature2);    
        }    
    
        OGRFeature::DestroyFeature(inputFeature1);    
        OGRGeometryFactory::destroyGeometry(cutGeometry);    
    }    
    
    // 关闭数据集    
    GDALClose(inputDataset1);    
    GDALClose(inputDataset2);    
    GDALClose(outputDataset);    
    
    // 关闭裁剪对话框    
    accept();    
    ui.InputlineEdit1->clear();    
    ui.InputlineEdit2->clear();    
    ui.OutputlineEdit->clear();    
}

5.3 矢量转栅格

功能点描述：将矢量图转化为栅格图层显示挂树，并输出tiff文件。

操作及截图：

当点击输入图层按钮时，弹出文件选择对话框，用户可以选择输入的Shapefile文件。选择后，将文件路径显示在界面上的InputEdit文本框中。然后使用GDAL库打开选定的Shapefile文件，并获取第一个图层的字段名称，将这些字段名称添加到界面上的下拉框(comboBox)中供用户选择。

图 93 获取输入图层

选择属性的combobox是获取输入图层的属性字段作为矢量转换字段，之后点击输出输出图层按钮弹出文件保存对话框，用户可以选择输出的栅格文件（TIFF文件）的保存路径。选择后，将文件路径显示在界面上的OutputEdit文本框中。

图 94 选定属性字段，确定输出路径

最后点击确定按钮时，开始进行转换操作。首先，获取输入和输出文件的路径以及用户选择的属性名称，注册GDAL驱动并打开输入Shapefile文件，获取输入Shapefile文件的第一个图层，并获取图层的空间范围和投影信息并创建一个栅格驱动，并创建一个指定大小的栅格数据集，用于保存转换后的栅格数据。根据图层的空间范围和大小，设置栅格数据集的地理转换参数（地理坐标系、像元大小等）。如果图层具有投影信息，将投影信息设置到栅格数据集中，设置好后创建一个包含要转换图层的数组，并设置转换选项。使用GDAL的GDALRasterizeLayers函数将图层转换为栅格数据。最后打开转换后的栅格文件，并将其添加到QGIS项目中进行显示。

图 95 矢量转栅格成果

核心代码：

void shptranstiff::on_InputpushButton_clicked() {  
    QString filePath = QFileDialog::getOpenFileName(nullptr, "Select shapefile", "", "Shapefiles (*.shp)");  
    ui.InputEdit->setText(filePath);  
    GDALAllRegister();  // 注册GDAL驱动  
    // 打开输入数据集  
    GDALDataset* inputDataset = (GDALDataset*)GDALOpenEx(filePath.toUtf8().constData(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr);  
    if (inputDataset == nullptr) {  
        QMessageBox::critical(this, "错误", "无法打开输入文件。");  
        return;  
    }  
  
    // 获取图层数量  
    int layerCount = inputDataset->GetLayerCount();  
  
    // 清空下拉框  
    ui.comboBox->clear();  
  
    // 获取选定图层索引  
    int selectedLayerIndex = 0; // 假设选择第一个图层  
    OGRLayer* selectedLayer = inputDataset->GetLayer(selectedLayerIndex);  
  
    // 获取字段数量  
    int fieldCount = selectedLayer->GetLayerDefn()->GetFieldCount();  
  
    // 用属性名称填充下拉框  
    for (int i = 0; i < fieldCount; i++) {  
        OGRFieldDefn* fieldDefn = selectedLayer->GetLayerDefn()->GetFieldDefn(i);  
        const char* attributeName = fieldDefn->GetNameRef();  
        ui.comboBox->addItem(QString::fromUtf8(attributeName));  
    }  
  
    // 关闭Shapefile  
    GDALClose(inputDataset);  
}  
  
void shptranstiff::on_OutputpushButton_clicked() {  
    QString saveFilePath = QFileDialog::getSaveFileName(nullptr, "Save tiff file", "", "files (*.tiff)");  
    ui.OutputEdit->setText(saveFilePath);  
}  
  
void shptranstiff::on_pushButton_clicked() {  
    QString inputPath = ui.InputEdit->text();  
    QString outputPath = ui.OutputEdit->text();  
    // 获取选定的属性名称  
    QString attributeName = ui.comboBox->currentText();  
  
    // 检查是否选择了属性名称  
    if (attributeName.isEmpty()) {  
        QMessageBox::critical(this, "错误", "请选择属性名称。");  
        return;  
    }  
    // 注册所有的驱动程序  
    GDALAllRegister();  
  
    // 打开矢量文件  
    GDALDataset* shpData = (GDALDataset*)GDALOpenEx(inputPath.toUtf8().constData(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr);  
    if (shpData == nullptr) {  
        QMessageBox::critical(this, "Error", "Failed to open vector file.");  
        return;  
    }  
  
    // 获取第一个图层  
    OGRLayer* shpLayer = shpData->GetLayer(0);  
    OGREnvelope env;  
    shpLayer->GetExtent(&env);//获取图层的坐标范围到env指向的内存中  
    int m_nImageWidth = 1440;  
    int m_nImageHeight = 720;  
    OGRSpatialReference* pOgrSRS = shpLayer->GetSpatialRef();//数据投影信息  
    char* pPrj = NULL;  
    if (pOgrSRS == NULL)  
    {  
        QMessageBox::critical(this, "Error", "没有投影信息");  
        m_nImageHeight = (int)env.MaxX;  
        m_nImageWidth = (int)env.MaxY;  
    }  
    else  
    {  
        pOgrSRS->exportToWkt(&pPrj);  
    }  
  
    // 创建栅格驱动  
    GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");  
    GDALDataset* poNewDS = poDriver->Create(outputPath.toStdString().c_str(), m_nImageWidth, m_nImageHeight, 1, GDT_Float32, NULL);  
    double adfGeoTransform[6];  
    adfGeoTransform[0] = env.MinX;//左上角经度  
    adfGeoTransform[1] = (env.MaxX - env.MinX) / m_nImageWidth;//像元宽度(影像在宽度上的分辨率)  
    adfGeoTransform[2] = 0; //如果影像是指北的, 这个参数的值为0  
    adfGeoTransform[3] = env.MaxY;//左上角纬度  
    adfGeoTransform[4] = 0; //如果影像是指北的, 这个参数的值为0。  
    adfGeoTransform[5] = (env.MinY - env.MaxY) / m_nImageHeight; //像元高度(影像在高度上的分辨率)  
    GDALSetGeoTransform(poNewDS, adfGeoTransform);  
    if (pOgrSRS != NULL)  
    {  
        poNewDS->SetProjection(pPrj);  
    }  
  
    int* pnbandlist = new int[1];  
    pnbandlist[0] = 1;  
  
    OGRLayerH* player;  
    player = new OGRLayerH[1];  
    player[0] = (OGRLayerH)shpLayer;  
  
    char** papszOptions = NULL;  
    papszOptions = CSLSetNameValue(papszOptions, "CHUNKSIZE", "1");  
    papszOptions = CSLSetNameValue(papszOptions, "ATTRIBUTE", attributeName.toUtf8().constData());  
  
    CPLErr err = GDALRasterizeLayers((GDALDatasetH)poNewDS, 1, pnbandlist, 1, player,NULL, NULL, NULL, papszOptions, NULL, NULL);  
  
    GDALClose(shpData);  
    GDALClose(poNewDS);  
    GDALDestroyDriverManager();  
    delete[]player;  
    delete[]pnbandlist;  
    close();  
    QMessageBox::information(this, "成功", "栅格文件转换完成。");  
  
    // 打开转换后的栅格文件  
    QFileInfo fi(outputPath);  
    if (!fi.exists()) { return; }  
    QString layerBaseName = fi.baseName(); // 图层名称  
    QgsRasterLayer* rasterLayer = new QgsRasterLayer(outputPath, layerBaseName);  
    if (!rasterLayer->isValid())  
    {  
        QMessageBox::critical(this, "Error", "Failed to load raster layer.");  
        return;  
    }  
  
    // 展示图层  
    QgsProject::instance()->addMapLayer(rasterLayer);  
      
}

5.4 栅格计算器

功能点描述：编写一个具有加、减、乘、除、对数、指数计算功能的栅格图层计算工具，计算结果以tiff格式的栅格文件输出。

操作及截图：

通过单击输入栅格图层按钮，调用on_InputpushButton_clicked()函数。函数中使用QFileDialog::getOpenFileName()函数打开文件对话框，让用户选择输入的TIFF文件。选择完成后，将文件路径显示在inputedit文本框中。然后使用GDAL库打开输入的TIFF图像文件，并将其作为GDAL数据集(GDALDataset)存储在inputDataset指针中。

之后通过组合框comboBox选择要计算的波段号。波段号从1开始，通过循环将波段号添加到组合框中。在计算表达式中输入进行计算的表达式；选择保存计算结果的位置：通过单击"OutputpushButton"按钮，调用on_OutputpushButton_clicked()函数。函数中使用QFileDialog::getSaveFileName()函数打开文件对话框，让用户选择保存计算结果的位置和文件名。选择完成后，将文件路径显示在outputedit文本框中。

最后点击"计算"按钮：通过单击"pushButton"按钮，调用on_pushButton_clicked()函数。函数中获取输入文件路径、输出文件路径和计算表达式。然后使用GDAL库打开输入的TIFF图像文件，并创建输出的TIFF图像文件。

图 96 栅格计算器演示界面

图 97 进行计算的源图层

图 98 栅格计算后的成果

核心代码：

void tiffcalculator::on_InputpushButton_clicked() {  
    QString filePath = QFileDialog::getOpenFileName(nullptr, "Select tiff", "", "Shapefiles (*.img *.tif *.tiff)");  
    ui.inputedit->setText(filePath);  
    GDALAllRegister();  // 注册GDAL驱动  
    // 打开输入数据集  
    GDALDataset* inputDataset = (GDALDataset*)GDALOpen(filePath.toUtf8().constData(), GA_ReadOnly);  
    if (inputDataset == nullptr) {  
        QMessageBox::critical(this, "错误", "无法打开输入文件。");  
        return;  
    }  
  
    // 使用波段号填充组合框  
    int bandCount = inputDataset->GetRasterCount();  
    ui.comboBox->clear();  
    for (int i = 1; i <= bandCount; ++i) {  
        ui.comboBox->addItem(QString::number(i));  
    }  
  
    GDALClose(inputDataset);  
}  
  
void tiffcalculator::on_OutputpushButton_clicked() {  
    QString saveFilePath = QFileDialog::getSaveFileName(nullptr, "Save tiff file", "", "files (*.tiff)");  
    ui.outputedit->setText(saveFilePath);  
}  
  
  
void tiffcalculator::on_pushButton_clicked() {  
    QString inputPath = ui.inputedit->text();  
    QString outputPath = ui.outputedit->text();  
    QString expression = ui.calculatelineEdit->text();  
    int selectedBand = ui.comboBox->currentText().toInt(); // 获取选择的波段号  
    GDALAllRegister();  // 注册GDAL驱动  
  
    // 打开输入图层  
    GDALDataset* inputDataset = (GDALDataset*)GDALOpen(inputPath.toUtf8().constData(), GA_ReadOnly);  
    if (inputDataset == nullptr) {  
        QMessageBox::critical(this, "错误", "打开输入图层失败");  
        return;  
    }  
  
    // 创建输出图层  
    GDALDriver* outputDriver = GetGDALDriverManager()->GetDriverByName("GTiff");  
    GDALDataset* outputDataset = outputDriver->Create(outputPath.toUtf8().constData(),inputDataset->GetRasterXSize(),inputDataset->GetRasterYSize(),1,GDT_Float32,nullptr);  
    if (outputDataset == nullptr) {  
        QMessageBox::critical(this, "错误", "创建输出图层失败");  
        GDALClose(inputDataset);  
        return;  
    }  
  
    // 获取输入图层的宽度和高度  
    int width = inputDataset->GetRasterXSize();  
    int height = inputDataset->GetRasterYSize();  
  
    // 读取输入图层数据  
    float* inputData = new float[width * height];  
    inputDataset->GetRasterBand(selectedBand)->RasterIO(GF_Read, 0, 0, width, height, inputData, width, height, GDT_Float32, 0, 0);  
      
    // 计算  
    float* outputData = new float[width * height];  
    for (int i = 0; i < width * height; ++i) {  
        float x = inputData[i];  
        float result = evaluateExpression(expression, x);  
        if (std::isnan(result)) {  
            QMessageBox::critical(this, "错误", "无效的表达式");  
            delete[] inputData;  
            delete[] outputData;  
            GDALClose(inputDataset);  
            GDALClose(outputDataset);  
            return;  
        }  
        outputData[i] = result;  
    }  
  
    // 将计算结果写入输出图层  
    outputDataset->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, width, height, outputData, width, height, GDT_Float32, 0, 0);  
  
    // 释放内存并关闭文件  
    delete[] inputData;  
    delete[] outputData;  
    GDALClose(inputDataset);  
    GDALClose(outputDataset);  
    close();  
    QMessageBox::information(this, "成功", "计算完成！");  
    // 打开计算后的栅格文件  
    QFileInfo fi(outputPath);  
    if (!fi.exists()) { return; }  
    QString layerBaseName = fi.baseName(); // 图层名称  
    QgsRasterLayer* rasterLayer = new QgsRasterLayer(outputPath, layerBaseName);  
    if (!rasterLayer->isValid())  
    {  
        QMessageBox::critical(this, "错误", "打开文件失败");  
        return;  
    }  
  
    // 展示图层  
    QgsProject::instance()->addMapLayer(rasterLayer);  
}