gdal test

https://blog.csdn.net/hb_programmer/article/details/81807699

 

gdal/ogr是一个光栅和矢量地理空间数据格式的翻译库，由开源地理空间基金会在开源许可下发布。作为一个库，它为所有支持的格式的调用应用程序提供单个光栅抽象数据模型和单个矢量抽象数据模型。它还附带了各种有用的命令行实用程序，用于数据转换和处理。新闻页面描述了2018年9月的gdal/ogr 2.3.2版本。

1. 

#include "gdal.h"
#include "gdal_priv.h"
#include "iostream"
using namespace std;

int main()
{
    //注册文件格式
    GDALAllRegister();

    const char* pszFile = "4.tif";
    //使用只读方式打开图像
    GDALDatasetH poDataset = GDALOpen(pszFile, GA_ReadOnly);
    GDALDataset *poSrcDS = GDALDataset::FromHandle(poDataset);
    if (poDataset == NULL)
    {

        printf("File：%s不能打开", pszFile);
        return 0;
    }
    //输出图像的格式信息                                                                                                                                                                                                                                                                                                                  
    printf("Driver:%s/%s\n",poSrcDS->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME),
        poSrcDS->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));//poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME)
    //输出图像的大小和波段个数
    printf("Size is %dx%dx%d\n", poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(), poSrcDS->GetRasterCount());
    //输出图像的投影信息
    if (poSrcDS->GetProjectionRef() != NULL)
        printf("Projection is '%s'\n",poSrcDS->GetProjectionRef());
    //输出图像的坐标和分辨率信息
    double adfGeoTransform[6];
    if (poSrcDS->GetGeoTransform(adfGeoTransform)==CE_None)
    {
        printf("Origin=(%.6f,%.6f)\n", adfGeoTransform[0], adfGeoTransform[3]);
        printf("PixelSize=(%.6f,%.6f)\n", adfGeoTransform[1], adfGeoTransform[5]);
    }
    //读取第一个波段
    GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
    //获取该波段的最大值最小值，如果获取失败，则进行统计
    int bGotMin, bGotMax;
    double adfMinMax[2];
    adfMinMax[0] = poBand->GetMinimum(&bGotMin);
    adfMinMax[1] = poBand->GetMaximum(&bGotMax);
    if (!(bGotMin&&bGotMax))
        GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
    printf("Min=%.3fd,Max=%.3f\n", adfMinMax[0], adfMinMax[1]);
    int nXSize = poBand->GetXSize();
    float *pafScanline = new float[nXSize];
    //读取图像的第一行数据
    poBand->RasterIO(GF_Read, 0, 0, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
    delete[]pafScanline;
    //关闭文件
    GDALClose((GDALDatasetH)poDataset);
    cin.get();
    //return 0;
}

2.

//构造OGRSpatialReference对象
void ConstructOSR()
{
    //定义一个OGRSpatialReference对象
    OGRSpatialReference oSRS;
    //下面使用三种方式对该对象进行初始化，设置为WGS84的地理坐标系
    //第一种方式
    oSRS.SetGeogCS("My geographic coordinate system", "WGS_1984", "My WGS84 Spheroid", SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
        "Greenwich", 0.0, "degree", 0.0174532925199433);
    //第二种方式使用常用名称进行构造
    oSRS.SetWellKnownGeogCS("WGS84");
    //第三种方式使用EPSG代码进行构造
    oSRS.SetWellKnownGeogCS("EPSG:4326");
    //第四种方式使用ESR的Prj文件进行构造
    //oSRS.SetFromUserInput("ESRI::C:\\Program Files(x86)\\ArcGIS\\Desktop10.0\\Coordinate Systems\\Geographic Coordinate Systems\\World\\WGS 1984.prj");
    //导出WKT格式
    char *pszWKT = NULL;
    oSRS.exportToWkt(&pszWKT);
    printf("%s\n\n", pszWKT);
    //导出美观的WKT格式
    char *pszPrettyWKT = NULL;
    oSRS.exportToPrettyWkt(&pszPrettyWKT);
    printf("%s\n", pszPrettyWKT);
}

//坐标转换
void CoordinateTransformation()
{
    const char* pszXian80 = "+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=39500000 +y_0=0 +ellps=IAU76 +towgs84=34.65192983,-69.97976937,-69.52875538,-0.56104022,-1.34050334,1.9067841,-0.27446825 +units=m _no_defs";
    OGRSpatialReference oXian80, *poLatLong;
    //使用PROJ.4格式字符串构造一个西安三度带投影坐标系统
    oXian80.SetFromUserInput(pszXian80);
    //获取该投影坐标系统中的地理坐标系统
    poLatLong = oXian80.CloneGeogCS();
    //构造一个从UTM投影坐标系统到地理坐标系统的转换关系
    OGRCoordinateTransformation *poTransform;
    poTransform = OGRCreateCoordinateTransformation(&oXian80, poLatLong);
    if (poTransform == NULL)
    {
        printf("创建坐标转换关系失败！\n");
        return;
    }
    double dX[2] = {39464667.861, 39458907.868};
    double dY[2] = {4441766.356, 444406.349};
    printf("转换前：\t1:(%f,%f)\n\t\t2:(%f,%f)\n", dX[0], dY[0], dX[1], dY[1]);
    if (!poTransform->Transform(2,dX,dY,NULL))
    {
        printf("坐标转换失败!\n");
        return;
    }
    printf("转换后:\t1:(%f,%f)\n\t\t2:(%f,%f)\n", dX[0], dY[0], dX[1], dY[1]);
}

3.

#include "ogrsf_frmts.h"
#include "iostream"
using namespace std;

bool LayerMethod(const char* pszSrcShp, const char* pszMethodShp, const char* pszDstShp, int iType)
{
  CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
  OGRRegisterAll();

  //打开数据
  //OGRDataSource *poSrcDS = OGRSFDriverRegistrar::GetOpenDS(pszSrcShp, FALSE);
  OGRDataSource poSrcDS;
  = OGRSFDriver::Open(pszSrcShp, FALSE);

}
void main()
{
  const char* pszSrcFile = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\城市绿地.shp";
  const char* pszMethodFile = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\111.shp";
  const char* pszOutFile0 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\0.shp";
  const char* pszOutFile1 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\1.shp";
  const char* pszOutFile2 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\2.shp";
  const char* pszOutFile3 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\3.shp";
  const char* pszOutFile4 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\4.shp";
  const char* pszOutFile5 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\5.shp";
  const char* pszOutFile6 = "F:\\本科毕业论文遥感GIS-滁州-植被适宜度评价\\666\\数据处理\\6.shp";

  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile0, 0);
  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile1, 1);
  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile2, 2);
  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile3, 3);
  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile4, 4);
  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile5, 5);
  LayerMethod(pszSrcFile, pszMethodFile, pszOutFile6, 6);
  cin.get();
}

3_2 test_ogr_shape.cpp

///////////////////////////////////////////////////////////////////////////////
//
// Project:  C++ Test Suite for GDAL/OGR
// Purpose:  Shapefile driver testing. Ported from ogr/ogr_shape.py.
// Author:   Mateusz Loskot <mateusz@loskot.net>
//
///////////////////////////////////////////////////////////////////////////////
// Copyright (c) 2006, Mateusz Loskot <mateusz@loskot.net>
// Copyright (c) 2010, Even Rouault <even dot rouault at mines-paris dot org>
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Library General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
// Boston, MA 02111-1307, USA.
///////////////////////////////////////////////////////////////////////////////

#include "gdal_unit_test.h"

#include "ogr_api.h"
#include "ogrsf_frmts.h"

#include <algorithm>
#include <iterator>
#include <string>
#include <vector>

namespace tut
{

    // Test data
    struct test_shape_data
    {
        OGRSFDriverH drv_;
        std::string drv_name_;
        std::string data_;
        std::string data_tmp_;

        test_shape_data()
            : drv_(nullptr), drv_name_("ESRI Shapefile")
        {
            drv_ = OGRGetDriverByName(drv_name_.c_str());

            // Compose data path for test group
            data_ = tut::common::data_basedir;
            data_tmp_ = tut::common::tmp_basedir;
        }
    };

    // Register test group
    typedef test_group<test_shape_data> group;
    typedef group::object object;
    group test_shape_group("OGR::Shape");

    // Test driver availability
    template<>
    template<>
    void object::test<1>()
    {
        ensure("OGR::Shape driver not available", nullptr != drv_);
    }

    // Test Create/Destroy empty directory datasource
    template<>
    template<>
    void object::test<2>()
    {
        // Try to remove tmp and ignore error code
        OGR_Dr_DeleteDataSource(drv_, data_tmp_.c_str());

        OGRDataSourceH ds = nullptr;
        ds = OGR_Dr_CreateDataSource(drv_, data_tmp_.c_str(), nullptr);
        ensure("OGR_Dr_CreateDataSource return NULL", nullptr != ds);

        OGR_DS_Destroy(ds);
    }

    // Create table from ogr/poly.shp
    template<>
    template<>
    void object::test<3>()
    {
        OGRErr err = OGRERR_NONE;

        OGRDataSourceH ds = nullptr;
        ds = OGR_Dr_CreateDataSource(drv_, data_tmp_.c_str(), nullptr);
        ensure("Can't open or create data source", nullptr != ds);

        // Create memory Layer
        OGRLayerH lyr = nullptr;
        lyr = OGR_DS_CreateLayer(ds, "tpoly", nullptr, wkbPolygon, nullptr);
        ensure("Can't create layer", nullptr != lyr);

        // Create schema
        OGRFieldDefnH fld = nullptr;

        fld = OGR_Fld_Create("AREA", OFTReal);
        err = OGR_L_CreateField(lyr, fld, true);
        OGR_Fld_Destroy(fld);
        ensure_equals("Can't create field", OGRERR_NONE, err);

        fld = OGR_Fld_Create("EAS_ID", OFTInteger);
        err = OGR_L_CreateField(lyr, fld, true);
        OGR_Fld_Destroy(fld);
        ensure_equals("Can't create field", OGRERR_NONE, err);

        fld = OGR_Fld_Create("PRFEDEA", OFTString);
        err = OGR_L_CreateField(lyr, fld, true);
        OGR_Fld_Destroy(fld);
        ensure_equals("Can't create field", OGRERR_NONE, err);

        // Check schema
        OGRFeatureDefnH featDefn = OGR_L_GetLayerDefn(lyr);
        ensure("Layer schema is NULL", nullptr != featDefn);
        ensure_equals("Fields creation failed", 3, OGR_FD_GetFieldCount(featDefn));

        // Copy ogr/poly.shp to temporary layer
        OGRFeatureH featDst = OGR_F_Create(featDefn);
        ensure("Can't create empty feature", nullptr != featDst);

        std::string source(data_);
        source += SEP;
        source += "poly.shp";
        OGRDataSourceH dsSrc = OGR_Dr_Open(drv_, source.c_str(), false);
        ensure("Can't open source layer", nullptr != dsSrc);

        OGRLayerH lyrSrc = OGR_DS_GetLayer(dsSrc, 0);
        ensure("Can't get source layer", nullptr != lyrSrc);

        OGRFeatureH featSrc = nullptr;
        while (nullptr != (featSrc = OGR_L_GetNextFeature(lyrSrc)))
        {
            err = OGR_F_SetFrom(featDst, featSrc, true);
            ensure_equals("Cannot set feature from source", OGRERR_NONE, err);

            err = OGR_L_CreateFeature(lyr, featDst);
            ensure_equals("Can't write feature to layer", OGRERR_NONE, err);

            OGR_F_Destroy(featSrc);
        }

        // Release and close resources
        OGR_F_Destroy(featDst);
        OGR_DS_Destroy(dsSrc);
        OGR_DS_Destroy(ds);
    }

    // Test attributes written to new table
    template<>
    template<>
    void object::test<4>()
    {
        OGRErr err = OGRERR_NONE;
        const int size = 5;
        const int expect[size] = { 168, 169, 166, 158, 165 };

        std::string source(data_tmp_);
        source += SEP;
        source += "tpoly.shp";
        OGRDataSourceH ds = OGR_Dr_Open(drv_, source.c_str(), false);
        ensure("Can't open layer", nullptr != ds);

        OGRLayerH lyr = OGR_DS_GetLayer(ds, 0);
        ensure("Can't get source layer", nullptr != lyr);

        err = OGR_L_SetAttributeFilter(lyr, "eas_id < 170");
        ensure_equals("Can't set attribute filter", OGRERR_NONE, err);

        // Prepare tester collection
        std::vector<int> list;
        std::copy(expect, expect + size, std::back_inserter(list));

        ensure_equal_attributes(lyr, "eas_id", list);

        OGR_DS_Destroy(ds);
    }

    // Test geometries written to new shapefile
    template<>
    template<>
    void object::test<5>()
    {
        // Original shapefile
        std::string orig(data_);
        orig += SEP;
        orig += "poly.shp";
        OGRDataSourceH dsOrig = OGR_Dr_Open(drv_, orig.c_str(), false);
        ensure("Can't open layer", nullptr != dsOrig);

        OGRLayerH lyrOrig = OGR_DS_GetLayer(dsOrig, 0);
        ensure("Can't get layer", nullptr != lyrOrig);

        // Copied shapefile
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH dsTmp = OGR_Dr_Open(drv_, tmp.c_str(), false);
        ensure("Can't open layer", nullptr != dsTmp);

        OGRLayerH lyrTmp = OGR_DS_GetLayer(dsTmp, 0);
        ensure("Can't get layer", nullptr != lyrTmp);

        // Iterate through features and compare geometries
        OGRFeatureH featOrig = OGR_L_GetNextFeature(lyrOrig);
        OGRFeatureH featTmp = OGR_L_GetNextFeature(lyrTmp);

        while (nullptr != featOrig && nullptr != featTmp)
        {
            OGRGeometryH lhs = OGR_F_GetGeometryRef(featOrig);
            OGRGeometryH rhs = OGR_F_GetGeometryRef(featTmp);

            ensure_equal_geometries(lhs, rhs, 0.000000001);

            // TODO: add ensure_equal_attributes()

            OGR_F_Destroy(featOrig);
            OGR_F_Destroy(featTmp);

            // Move to next feature
            featOrig = OGR_L_GetNextFeature(lyrOrig);
            featTmp = OGR_L_GetNextFeature(lyrTmp);
        }

        OGR_DS_Destroy(dsOrig);
        OGR_DS_Destroy(dsTmp);
    }

    // Write a feature without a geometry
    template<>
    template<>
    void object::test<6>()
    {
        // Create feature without geometry
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH ds = OGR_Dr_Open(drv_, tmp.c_str(), true);
        ensure("Can't open layer", nullptr != ds);

        OGRLayerH lyr = OGR_DS_GetLayer(ds, 0);
        ensure("Can't get layer", nullptr != lyr);

        OGRFeatureDefnH featDefn = OGR_L_GetLayerDefn(lyr);
        ensure("Layer schema is NULL", nullptr != featDefn);

        OGRFeatureH featNonSpatial = OGR_F_Create(featDefn);
        ensure("Can't create non-spatial feature", nullptr != featNonSpatial);

        int fldIndex = OGR_FD_GetFieldIndex(featDefn, "PRFEDEA");
        ensure("Can't find field 'PRFEDEA'", fldIndex >= 0);

        OGR_F_SetFieldString(featNonSpatial, fldIndex, "nulled");

        OGRErr err = OGR_L_CreateFeature(lyr, featNonSpatial);
        ensure_equals("Can't write non-spatial feature to layer", OGRERR_NONE, err);

        OGR_F_Destroy(featNonSpatial);
        OGR_DS_Destroy(ds);
    }

    // Read back the non-spatial feature and get the geometry
    template<>
    template<>
    void object::test<7>()
    {
        OGRErr err = OGRERR_NONE;

        // Read feature without geometry
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH ds = OGR_Dr_Open(drv_, tmp.c_str(), false);
        ensure("Can't open layer", nullptr != ds);

        OGRLayerH lyr = OGR_DS_GetLayer(ds, 0);
        ensure("Can't get layer", nullptr != lyr);

        err = OGR_L_SetAttributeFilter(lyr, "PRFEDEA = 'nulled'");
        ensure_equals("Can't set attribute filter", OGRERR_NONE, err);

        // Fetch feature without geometry
        OGRFeatureH featNonSpatial = OGR_L_GetNextFeature(lyr);
        ensure("Didn't get feature with null geometry back", nullptr != featNonSpatial);

        // Null geometry is expected
        OGRGeometryH nonGeom = OGR_F_GetGeometryRef(featNonSpatial);
        ensure("Didn't get null geometry as expected", nullptr == nonGeom);

        OGR_F_Destroy(featNonSpatial);
        OGR_DS_Destroy(ds);
    }

    // Test ExecuteSQL() results layers without geometry
    template<>
    template<>
    void object::test<8>()
    {
        const int size = 11;
        const int expect[size] = { 179, 173, 172, 171, 170, 169, 168, 166, 165, 158, 0 };

        // Open directory as a datasource
        OGRDataSourceH ds = OGR_Dr_Open(drv_, data_tmp_ .c_str(), false);
        ensure("Can't open datasource", nullptr != ds);

        std::string sql("select distinct eas_id from tpoly order by eas_id desc");
        OGRLayerH lyr = OGR_DS_ExecuteSQL(ds, sql.c_str(), nullptr, nullptr);
        ensure("Can't create layer from query", nullptr != lyr);

        // Prepare tester collection
        std::vector<int> list;
        std::copy(expect, expect + size, std::back_inserter(list));

        ensure_equal_attributes(lyr, "eas_id", list);

        OGR_DS_ReleaseResultSet(ds, lyr);
        OGR_DS_Destroy(ds);
    }

    // Test ExecuteSQL() results layers with geometry
    template<>
    template<>
    void object::test<9>()
    {
        // Open directory as a datasource
        OGRDataSourceH ds = OGR_Dr_Open(drv_, data_tmp_ .c_str(), false);
        ensure("Can't open datasource", nullptr != ds);

        std::string sql("select * from tpoly where prfedea = '35043413'");
        OGRLayerH lyr = OGR_DS_ExecuteSQL(ds, sql.c_str(), nullptr, nullptr);
        ensure("Can't create layer from query", nullptr != lyr);

        // Prepare tester collection
        std::vector<std::string> list;
        list.push_back("35043413");

        // Test attributes
        ensure_equal_attributes(lyr, "prfedea", list);

        // Test geometry
        const char* wkt = "POLYGON ((479750.688 4764702.000,479658.594 4764670.000,"
            "479640.094 4764721.000,479735.906 4764752.000,"
            "479750.688 4764702.000))";

        OGRGeometryH testGeom = nullptr;
        OGRErr err = OGR_G_CreateFromWkt((char**) &wkt, nullptr, &testGeom);
        ensure_equals("Can't create geometry from WKT", OGRERR_NONE, err);

        OGR_L_ResetReading(lyr);
        OGRFeatureH feat = OGR_L_GetNextFeature(lyr);
        ensure("Cannot fetch feature", nullptr != feat);

        ensure_equal_geometries(OGR_F_GetGeometryRef(feat), testGeom, 0.001);

        OGR_F_Destroy(feat);
        OGR_G_DestroyGeometry(testGeom);
        OGR_DS_ReleaseResultSet(ds, lyr);
        OGR_DS_Destroy(ds);
    }

    // Test spatial filtering
    template<>
    template<>
    void object::test<10>()
    {
        OGRErr err = OGRERR_NONE;

        // Read feature without geometry
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH ds = OGR_Dr_Open(drv_, tmp.c_str(), false);
        ensure("Can't open layer", nullptr != ds);

        OGRLayerH lyr = OGR_DS_GetLayer(ds, 0);
        ensure("Can't get layer", nullptr != lyr);

        // Set empty filter for attributes
        err = OGR_L_SetAttributeFilter(lyr, nullptr);
        ensure_equals("Can't set attribute filter", OGRERR_NONE, err);

        // Set spatial filter
        const char* wkt = "LINESTRING(479505 4763195,480526 4762819)";
        OGRGeometryH filterGeom = nullptr;
        err = OGR_G_CreateFromWkt((char**) &wkt, nullptr, &filterGeom);
        ensure_equals("Can't create geometry from WKT", OGRERR_NONE, err);

        OGR_L_SetSpatialFilter(lyr, filterGeom);

        // Prepare tester collection
        std::vector<int> list;
        list.push_back(158);

        // Test attributes
        ensure_equal_attributes(lyr, "eas_id", list);

        OGR_G_DestroyGeometry(filterGeom);
        OGR_DS_Destroy(ds);
    }

} // namespace tut

 

https://gdal.org/api/index.html