GDAL使用(二)
本篇说一下CreateCopy()和Create()
CreateCopy()是从一个源数据拷贝到一个新的数据中。
有关说明写到了代码的注释中。C++代码和Python代码实现的功能是一致的。
C++代码:
#include <iostream>
#include "gdal_priv.h"
#include "cpl_conv.h"
#include <cpl_string.h>
int main()
{
clock_t start = clock();
GDALAllRegister();
const char* srcDs = "E:\testdata\CZIOrigin.tiff";
const char* dstDs = "O:\landsatdata.tif";
const char* pszFormat = "GTiff";
//打开数据
GDALDataset* poSrcDs = (GDALDataset*)GDALOpen(srcDs, GA_ReadOnly);
//创建驱动
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if (poDriver == NULL)
exit(1);
//相关创建选项可以到栅格驱动程序中查看
char** papszOptions = NULL;
papszOptions = CSLSetNameValue(papszOptions, "TILED", "YES");
papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "LZW");
papszOptions = CSLSetNameValue(papszOptions, "NUM_THREADS", "ALL_CPUS");
//创建拷贝数据集
GDALDataset* poDstDs = poDriver->CreateCopy(dstDs, poSrcDs, FALSE, papszOptions, GDALTermProgress, NULL);
//释放资源
if (poDstDs != NULL)
{
GDALClose((GDALDatasetH)poDstDs);
}
GDALClose((GDALDatasetH)poSrcDs);
CSLDestroy(papszOptions);
clock_t end = clock();
double endtime = (double)(end - start) / CLOCKS_PER_SEC;
//以秒为单位输出时间
std::cout << "Total time:" << endtime << std::endl;
printf("over!");
system("pause");
}Python代码:
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 4 11:47:06 2021
@author: cyhu
"""
from osgeo import gdal
import time
time_start=time.time()
# 文件路径
src_file = "E:\testdata\landsatdata.tiff"
dst_file = "O:\landsatdata.tif"
# 打开数据
src_ds = gdal.Open(src_file,gdal.GA_ReadOnly)
if not src_ds:
print("open failed!")
driver = gdal.GetDriverByName("GTiff")
# 相关的创建选项可以到栅格驱动程序中查看
dst_ds = driver.CreateCopy(dst_file,src_ds,strict = 0,
options=["TILED=YES", "COMPRESS=LZW","NUM_THREADS=ALL_CPUS"])
src_ds = None
dst_ds = None
time_end = time.time()
print("time ",time_end - time_start)Create()本身并不复杂,但是需要进行一些元数据的写入,所以代码量稍多一些。
C++代码:
#include <iostream>
#include "gdal_priv.h"
#include "cpl_conv.h"
#include <cpl_string.h>
int main()
{
//clock_t start = clock();
GDALAllRegister();
//const char* srcDs = "E:\testdata\CZIOrigin.tiff";
const char* dstDs = "O:\landsatdata.tif";
const char* pszFormat = "GTiff";
//创建栅格驱动程序
GDALDriver* poDriver = (GDALDriver*)GDALGetDriverByName(pszFormat);
//创建栅格
char** papszOptions = NULL;
GDALDataset* poDstDS = poDriver->Create(dstDs, 512, 512, 1, GDT_Byte, papszOptions);
/* 创建上面的栅格数据后写入元数据和栅格数据 */
//要写入栅格的数据块
GByte abyRaster[512 * 512];
//创建仿射变换参数
double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
//创建OGR空间参考
OGRSpatialReference oSRS;
char* pszSRS_WKT = NULL;
//设置数据集的变换参数
poDstDS->SetGeoTransform(adfGeoTransform);
//设置投影分带
oSRS.SetUTM(11, TRUE);
//设置基准面
oSRS.SetWellKnownGeogCS("NAD27");
//输出到WKT格式
oSRS.exportToWkt(&pszSRS_WKT);
//设置投影
poDstDS->SetProjection(pszSRS_WKT);
//获取波段
GDALRasterBand* poBand = poDstDS->GetRasterBand(1);
//写入数据
poBand->RasterIO(GF_Write, 0, 0, 512, 512, abyRaster, 512, 512, GDT_Byte, 0, 0);
//释放资源
CPLFree(pszSRS_WKT);
GDALClose((GDALDatasetH)poDstDS);
std::cout << "over!" << std::endl;
system("pause");
return 0;
}Python代码:
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 4 17:21:33 2021
@author: cyhu
"""
from osgeo import gdal
from osgeo import osr
import numpy
# 获取栅格驱动
driver = gdal.GetDriverByName("GTiff")
dst_ds = "O:\test.tif"
# 创建栅格
dataset = driver.Create(dst_ds, xsize = 512, ysize = 512, bands = 1, eType = gdal.GDT_Byte)
# 设置仿射变换参数
dataset.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
# 设置空间参考
spatialref = osr.SpatialReference()
spatialref.SetUTM(11,1)
spatialref.SetWellKnownGeogCS("NAD27")
dataset.SetProjection(spatialref.ExportToWkt())
# 创建要写入的数据
raster = numpy.zeros((512, 512), dtype=numpy.uint8)
# 写入数据
dataset.GetRasterBand(1).WriteArray(raster)
# 释放资源
dataset = None


鲁公网安备 37148202000241号