C++使用GDAL输出遥感影像的.shp文件

C++使用GDAL输出遥感影像的.shp文件

注意：本文使用的GDAL版本为2.0.1

工作中遇到了一个小问题，需要输出一个指定影像的shp文件和prj文件。其中shp的数据格式是指定好的。

经过一番搜索，总算找到了可用的方法，经过修改，最终达到了我的要求。这里把最终版本的代码分享出来，供大家参考。

---

1. 从Tiff中获取投影信息

首先，使用GDAL读取影像文件，获取其中的投影信息。当然，需要文件中已经包含有这些信息，否则读取到的字符串为NULL。

另外，GDAL也提供了获取图像地理转换六参数的方法，通过GetGeoTransform函数将转换参数读取到浮点型指针中，指针记得提前分配好内存，否则会报未初始化的错误。

#include <iostream>
#include <string>
// For GDAL
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
using namespace std;
int main (int argc, char* argv[])
{
    // 注册所需驱动
    GDALAllRegister();
    string proj = "", filepath = "F:/1.tiff";
    double* geotrans = new double[6];
    GDALDataset * poDataset = (GDALDataset *) GDALOpen (filepath.c_str(), GA_ReadOnly);
    
    // 打开影像不能为空
    if (poDataset == NULL)
    {
        cout << "Read tiff error!";
        return 1;
    }
    
    // 获取一个字符串，即使我们所需的影像投影信息
    proj = poDataset->GetProjectionRef();
    // 获取地理转换六参数
    poDataset->GetGeoTransform (geotrans);
    GDALClose (poDataset);
    return 0;
}

上述代码获取到的投影信息如下：

GEOGCS[“GCS_WGS_1984”,DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137,298.2572235630016]],PRIMEM[“Greenwich”,0],UNIT[“Degree”,0.017453292519943295]]

---

2. Tiff输出对应的.shp.shx.tfw文件

直接上代码，文件输出到shpfilepath路径下。

#include <iostream>
#include <string>
// For GDAL
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
using namespace std;

int main (int argc, char* argv[])
{
    // 注册所需驱动
    GDALAllRegister();
    OGRRegisterAll();
    const char *pszDriverName = "ESRI Shapefile";
    string shpfilepath = "F:\\1.shp";
    GDALDriver *poDriver;
    GDALAllRegister();
    poDriver = GetGDALDriverManager()->GetDriverByName (pszDriverName);
    
    if (poDriver == NULL)
    {
        cout << "Driver not available." << endl;
        return 1;
    }
    
    GDALDataset *poDS;
    poDS = poDriver->Create (shpfilepath.c_str(), 0, 0, 0, GDT_Unknown, NULL);
    
    if (poDS == NULL)
    {
        cout << "Creation of output file failed." << endl;;
        return 1;
    }
    
    OGRLayer *poLayer;
    poLayer = poDS->CreateLayer ("Polygon", NULL, wkbPolygon, NULL);
    
    if (poLayer == NULL)
    {
        cout <<  "Layer creation failed." << endl;
        return 1;
    }
    
    double TLlat = 0, TLlon = 0;// 影像左上角点经纬度
    double TRlat = 0, TRlon = 0.1;// 影像右上角点经纬度
    double BLlat = 0.1, BLlon = 0;// 影像左下角点经纬度
    double BRlat = 0.1, BRlon = 0.1;// 影像右下角点经纬度
    OGRFieldDefn f_Index ("index", OFTInteger);// 添加整型字段index
    OGRFieldDefn f_NAME ("NAME", OFTString);// 添加字符串字段NAME
    OGRFieldDefn f_TopLeftLat ("TopLeftLat", OFTReal);// 添加浮点型字段TopLeftLat
    OGRFieldDefn f_TopLeftLon ("TopLeftLon", OFTReal);
    OGRFieldDefn f_TopRightLa ("TopRightLa", OFTReal);
    OGRFieldDefn f_TopRightLo ("TopRightLo", OFTReal);
    OGRFieldDefn f_BotLeftLat ("BotLeftLat", OFTReal);
    OGRFieldDefn f_BotLeftLon ("BotLeftLon", OFTReal);
    OGRFieldDefn f_BotRightLa ("BotRightLa", OFTReal);
    OGRFieldDefn f_BotRightLo ("BotRightLo", OFTReal);
    // 设置字段宽度
    f_Index.SetWidth (32);
    f_NAME.SetWidth (128);
    f_TopLeftLat.SetWidth (23);
    f_TopLeftLat.SetPrecision (15);// 对于浮点型字段，需要设置其精度；默认为0时，只存储整数部分。
    f_TopLeftLon.SetWidth (23);
    f_TopLeftLon.SetPrecision (15);
    f_TopRightLa.SetWidth (23);
    f_TopRightLa.SetPrecision (15);
    f_TopRightLo.SetWidth (23);
    f_TopRightLo.SetPrecision (15);
    f_BotLeftLat.SetWidth (23);
    f_BotLeftLat.SetPrecision (15);
    f_BotLeftLon.SetWidth (23);
    f_BotLeftLon.SetPrecision (15);
    f_BotRightLa.SetWidth (23);
    f_BotRightLa.SetPrecision (15);
    f_BotRightLo.SetWidth (23);
    f_BotRightLo.SetPrecision (15);
    // 添加字段
    poLayer->CreateField (&f_Index);
    poLayer->CreateField (&f_NAME);
    poLayer->CreateField (&f_TopLeftLat);
    poLayer->CreateField (&f_TopLeftLon);
    poLayer->CreateField (&f_TopRightLa);
    poLayer->CreateField (&f_TopRightLo);
    poLayer->CreateField (&f_BotLeftLat);
    poLayer->CreateField (&f_BotLeftLon);
    poLayer->CreateField (&f_BotRightLa);
    poLayer->CreateField (&f_BotRightLo);
    OGRFeature *poFeature;
    poFeature = OGRFeature::CreateFeature (poLayer->GetLayerDefn());
    poFeature->SetField ("index", 0);
    poFeature->SetField ("NAME", "ID");
    poFeature->SetField ("TopLeftLat", TLlat);
    poFeature->SetField ("TopLeftLon", TLlon);
    poFeature->SetField ("TopRightLa", TRlat);
    poFeature->SetField ("TopRightLo", TRlon);
    poFeature->SetField ("BotLeftLat", BLlat);
    poFeature->SetField ("BotLeftLon", BLlon);
    poFeature->SetField ("BotRightLa", BRlat);
    poFeature->SetField ("BotRightLo", BRlon);
    // 定义一个闭环，用来表示矩形图像范围
    // 这里要注意矩形的四个角点的添加顺序，要按顺时针或逆时针添加
    OGRLinearRing ring;
    ring.addPoint (TLlon, TLlat);
    ring.addPoint (TRlon, TRlat);
    ring.addPoint (BRlon, BRlat);
    ring.addPoint (BLlon, BLlat);
    ring.closeRings();
    OGRPolygon polygon;
    polygon.addRing (&ring);
    poFeature->SetGeometry (&polygon);

    // 添加Feature到图层中
    if (poLayer->CreateFeature (poFeature) != OGRERR_NONE)
    {
        cout <<  "Failed to create feature in shapefile." << endl;
        return 1;
    }

	// 关闭并存储图层
    OGRFeature::DestroyFeature (poFeature);
    GDALClose (poDS);
    return 0;
}

把输出的.shp文件用ArcMap打开，就可以看到刚刚定义的属性了。其中的字段信息也可以根据实际需要进行定义。

参考

1.结合C++和GDAL实现shapefile（shp)文件的创建和写入