使用GDAL/OGR读写矢量文件

感觉GIS中矢量相关内容还是挺庞杂的，并且由于版本迭代的关系，使用GDAL/OGR读写矢量的资料也有点不太一样。这里总结了一个读写矢量的示例，实现代码如下：

代码语言：javascript复制

#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文件：

gdal gis

0 人点赞