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