使用GDAL/OGR讀寫矢量文件

感受GIS中矢量相關內容仍是挺龐雜的，而且因爲版本迭代的關係，使用GDAL/OGR讀寫矢量的資料也有點不太同樣。這裏總結了一個讀寫矢量的示例，實現代碼以下：

#include <iostream>

#include <gdal/ogrsf_frmts.h>

using namespace std;

bool ReadDXF(string filePath, vector<vector<OGRPoint>>& vertexPoint)
{   
    GDALDataset *poDS = (GDALDataset*)GDALOpenEx(filePath.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (!poDS)
    {
        printf("沒法讀取該文件，試檢查格式是否正確！");
        return false;
    }
    if (poDS->GetLayerCount()<1)
    {
        printf("該文件的層數小於1，試檢查格式是否正確！");
        return false;
    }

    OGRLayer  *poLayer = poDS->GetLayer(0); //讀取層
    poLayer->ResetReading();
    
    OGRFeature *poFeature;
    while ((poFeature = poLayer->GetNextFeature()) != NULL)
    {
        OGRGeometry *pGeo = poFeature->GetGeometryRef();
        OGRwkbGeometryType pGeoType = pGeo->getGeometryType();

        if (pGeoType == wkbLineString || pGeoType == wkbLineString25D)
        {
            OGRLinearRing  *pCurve = (OGRLinearRing*)pGeo;
            if (pCurve->getNumPoints() < 1)
            {
                continue;
            }

            vector<OGRPoint> pl;
            for (int i = 0; i<pCurve->getNumPoints(); i++)
            {
                OGRPoint point;
                pCurve->getPoint(i, &point);                        
                pl.push_back(point);
            }
            vertexPoint.push_back(pl);
        }

        ////        
        //OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();        
        //int n = poFDefn->GetFieldCount(); //得到字段的數目，不包括前兩個字段（FID,Shape);
        //for (int iField = 0; iField <n; iField++)
        //{         
        //    //輸出每一個字段的值
        //    cout << poFeature->GetFieldAsString(iField) << "    ";            
        //}
        //cout << endl;   

        OGRFeature::DestroyFeature(poFeature);
    }

    GDALClose(poDS);
    poDS = nullptr;

    return true;
}

bool WriteShp(string filePath, vector<vector<OGRPoint>> vertexPoint)
{
    //建立
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");
    if (!driver)
    {
        printf("Get Driver ESRI Shapefile Error！\n");
        return false;
    }
    
    GDALDataset* dataset = driver->Create(filePath.c_str(), 0, 0, 0, GDT_Unknown, NULL);
    OGRLayer* poLayer = dataset->CreateLayer("houseType", NULL, wkbPolygon, NULL);
    
    //建立屬性字段
    {
        // 字符串
        OGRFieldDefn oField1("名稱", OFTString);
        oField1.SetWidth(8);
        if (poLayer->CreateField(&oField1) != OGRERR_NONE) {
            printf("Creating Name field failed.\n"); return FALSE;
        }

        // 浮點數
        OGRFieldDefn oField2("面積", OFTReal);
        oField2.SetPrecision(3);
        if (poLayer->CreateField(&oField2) != OGRERR_NONE) {
            printf("Creating Name field failed.\n"); return FALSE;
        }

        // 整型
        OGRFieldDefn oField3("結點數", OFTInteger);
        if (poLayer->CreateField(&oField3) != OGRERR_NONE) {
            printf("Creating Name field failed.\n"); return FALSE;
        }
    }

    //建立特徵
    for (auto& iter : vertexPoint)
    {
        OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn());
        
        OGRLinearRing ogrring;
        int pNum = (int)iter.size();
        ogrring.setNumPoints(pNum);
        for (int i = 0; i < iter.size(); i++)
        {
            ogrring.setPoint(i, iter[i].getX(), iter[i].getY(), iter[i].getZ());
            //cout << iter[i].x() << '\t' << iter[i].y() << '\t' << iter[i].z() << endl;
        }
        //cout << "-----------------------------\n";

        OGRPolygon polygon;
        polygon.addRing(&ogrring);
        poFeature->SetGeometry(&polygon);       

        poFeature->SetField("名稱", "多邊形");
        poFeature->SetField("面積", polygon.get_Area());
        poFeature->SetField("結點數", pNum);

        if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
        {
            printf("Failed to create feature in shapefile.\n");
            return false;
        }
    }
    
    //釋放
    GDALClose(dataset);
    dataset = nullptr;
    //GDALDestroyDriverManager();

    return true;
}

int main()
{
    GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路徑
    CPLSetConfigOption("SHAPE_ENCODING", "");  //解決中文亂碼問題
        
    string filePath = "D:/2.dxf";
    vector<vector<OGRPoint>> vertexPoint;
    if (!ReadDXF(filePath, vertexPoint))
    {
        return 1;
    }

    string newPath = "C:/Users/charlee/Desktop/SHP/dst.shp";
    WriteShp(newPath, vertexPoint);

    return 0;
}

在這個示例中，讀取一個DXF文件中的線（環）特徵，將其轉換成面，而後保存在一個SHP中。同時，還給該SHP文件寫入了相應的屬性字段。

讀取的DXF文件：

建立並保存的SHP文件：