当前位置：首页 > C++ > 正文内容

FY4B warp

admin2年前 (2024-01-03)C++16583

```cpp
#include <iostream>
#include <string>
#include <gdal_priv.h>
#include "gdalwarper.h"

using std::string;

void warpdata(GDALDatasetH hSrcDS, GDALDatasetH hDstDS)
{

// Open input and output files.
	GDALAllRegister();

// Setup warp options.
	GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
	psWarpOptions->hSrcDS = hSrcDS;
	psWarpOptions->hDstDS = hDstDS;
	psWarpOptions->nBandCount = 1;
	psWarpOptions->panSrcBands = (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
	psWarpOptions->panSrcBands[0] = 1;
	//psWarpOptions->panSrcBands[1] = 2;
	//psWarpOptions->panSrcBands[2] = 3;
	//psWarpOptions->panSrcBands[3] = 4;
	psWarpOptions->panDstBands = (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
	psWarpOptions->panDstBands[0] = 1;
	//psWarpOptions->panDstBands[1] = 2;
	//psWarpOptions->panDstBands[2] = 3;
	//psWarpOptions->panDstBands[3] = 4;
	psWarpOptions->pfnProgress = GDALTermProgress;

// Establish reprojection transformer.
	const char* wktName = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
	OGRSpatialReference* spatial = new OGRSpatialReference();
	spatial->importFromProj4("+proj=geos +no_defs +a=6378137.0 +b=6356752.3 +lon_0=133.0 +h=35785863");
	char* wkt = NULL;
	spatial->exportToWkt(&wkt);
	char** papszOptions = NULL;
	papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", wkt);
	papszOptions = CSLSetNameValue(papszOptions, "DST_SRS", wktName);
	char** papszWarpOptions = NULL;
	papszWarpOptions = CSLSetNameValue(papszWarpOptions, "NUM_THREADS", "ALL_CPUS");
	psWarpOptions->papszWarpOptions = papszWarpOptions;

//papszOptions = CSLSetNameValue(papszOptions, "SRC_METHOD", "NO_GEOTRANSFORM");
	psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer2(hSrcDS, hDstDS, papszOptions);
	psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

// Initialize and execute the warp operation.
	GDALWarpOperation oOperation;
	oOperation.Initialize(psWarpOptions);
	oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));
	GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
	GDALDestroyWarpOptions(psWarpOptions);
}

int main(int argc, char* argv[])
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");//设置支持中文路径
	//string FDIdataPath = argv[1];
	//string LDAdataPath = argv[2];

//// 打开HDF文件
	//GDALDataset* hDataset = (GDALDataset*)GDALOpen("F:\\FY-4\\FY4A\\FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20220815063000_20220815063417_1000M_V0001.HDF", GA_ReadOnly);
	//if (hDataset == nullptr) {
	//    std::cerr << "无法打开HDF文件" << std::endl;
	//    return -1;
	//}
	//// 获取HDF文件的元数据
	//const char* metadata = hDataset->GetMetadataItem("Begin_Line_Number");
	//if (metadata != nullptr) {
	//    std::cout << "HDF元数据: " << metadata << std::endl;
	//}
	//else {
	//    std::cout << "无法获取HDF元数据" << std::endl;
	//}

//
	//GDALClose(hDataset);

string FDIdataPath = "Z:\\PersonalFolder\\cyhu\\FY4B大气校正样例数据\\FY4B-_AGRI--_N_DISK_1330E_L1-_FDI-_MULT_NOM_20221011030000_20221011031459_1000M_V0001.HDF";

//double LonMin = 51.9101;// std::stod(argv[5]);
	//double LonMax = 214.2701;// std::stod(argv[6]);
	//double LatMin = -80.9922;// std::stod(argv[7]);
	//double LatMax = 81.0578;//std::stod(argv[8]);
	////计算图像的长和高，步长（分辨率）为0.01
	//double step = 0.01;

//int xsize = (int)((LonMax - LonMin) / step);
	//int ysize = (int)((LatMax - LatMin) / step);

string b1Path = "HDF5:\"" + FDIdataPath + "\"://Data/NOMChannel01";
	string b1CalibriPath = "HDF5:\"" + FDIdataPath + "\"://Calibration/CALChannel01";
	string b2Path = "HDF5:\"" + FDIdataPath + "\"://Data/NOMChannel02";
	string b2CalibriPath = "HDF5:\"" + FDIdataPath + "\"://Calibration/CALChannel02";
	string b3Path = "HDF5:\"" + FDIdataPath + "\"://Data/NOMChannel03";
	string b3CalibriPath = "HDF5:\"" + FDIdataPath + "\"://Calibration/CALChannel03";
	
	GDALDataset* dst1 = (GDALDataset*)GDALOpen(b1Path.c_str(), GA_ReadOnly);
	//GDALDataset* dst1Calibri = (GDALDataset*)GDALOpen(b1CalibriPath.c_str(), GA_ReadOnly);

//GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	//GDALDataset* outds1 = driver->Create("F:\\testimgs\\rasterfiles\\FY4BGEOREF.tif", xsize, ysize, 3, GDT_Int16, 0);
	const char* wktName = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
	//double geotransform[6] = { LonMin, 0.01, 0.0, LatMax, 0.0, -0.01 };
	//outds1->SetProjection(wktName);
	//outds1->SetGeoTransform(geotransform);

//int originXsize = dst1->GetRasterXSize();
	//int originYsize = dst1->GetRasterYSize();

//int calibriDataXsize = dst1Calibri->GetRasterXSize();
	//int calibriDataYsize = dst1Calibri->GetRasterYSize();
	//
	//unsigned short* origionData = new unsigned short[originXsize * originYsize];//用于保存原始数据，每次原始数据都是读取到这个数组里面
	//float* calibridata = new float[calibriDataXsize * calibriDataYsize]; //用于保存定标查找表数据，每次原始数据都是读取到这个数组里面
	//dst1->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, originXsize, originYsize, origionData, originXsize, originYsize, GDT_UInt16, 0, 0);
	//dst1Calibri->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, calibriDataXsize, calibriDataYsize, calibridata, calibriDataXsize, calibriDataYsize, GDT_Float32, 0, 0);
	//short* outb1data = new short[ysize * xsize];//b band
	//short* lutcol = new short[ysize * xsize];//lut
	//short* lutline = new short[ysize * xsize];//lut
	//for (int i = 0; i < ysize; i++)
	//{
	//	for (int j = 0; j < xsize; j++)
	//	{
	//		int line = 0;
	//		int col = 0;
	//		LatLonToRC(LonMin + step * j, LatMax - i * step, col, line);
	//		lutcol[i * xsize + j] = col;
	//		lutline[i * xsize + j] = line;
	//	}

//	if (i * 1.0 / ysize > 2 && i * 1.0 / ysize < 3)
	//	{
	//		std::cout << i * 1.0 / ysize << ".";
	//	}
	//	
	//}

//for (int i = 0; i < ysize; i++)
	//{
	//	for (int j = 0; j < xsize; j++)
	//	{
	//		short line = lutline[i * xsize + j];
	//		short col = lutcol[i * xsize + j];
	//		//LatLonToRC(LonMin + step * j, LatMax - i * step, col, line);
	//		//outb1data[i * xsize + j] = origionData[originXsize * line + col];
	//		int index = origionData[originXsize * line + col];
	//		if (index != 65535)
	//		{
	//			outb1data[i * xsize + j] = short(calibridata[index] * 1000);
	//		}
	//		else
	//		{
	//			outb1data[i * xsize + j] = 0;
	//		}
	//		
	//	}
	//}
	//outds1->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, xsize, ysize, outb1data, xsize, ysize, GDT_Int16, 0, 0, 0);

GDALDataset* dst2 = (GDALDataset*)GDALOpen(b2Path.c_str(), GA_ReadOnly);
	GDALDataset* dst2Calibri = (GDALDataset*)GDALOpen(b2CalibriPath.c_str(), GA_ReadOnly);
	//dst2->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, originXsize, originYsize, origionData, originXsize, originYsize, GDT_UInt16, 0, 0);
	//dst2Calibri->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, calibriDataXsize, calibriDataYsize, calibridata, calibriDataXsize, calibriDataYsize, GDT_Float32, 0, 0);
	//for (int i = 0; i < ysize; i++)
	//{
	//	for (int j = 0; j < xsize; j++)
	//	{
	//		int line = lutline[i * xsize + j];
	//		int col = lutcol[i * xsize + j];
	//		//LatLonToRC(LonMin + step * j, LatMax - i * step, col, line);
	//		//outb1data[i * xsize + j] = origionData[originXsize * line + col];
	//		int index = origionData[originXsize * line + col];
	//		outb1data[i * xsize + j] = short(calibridata[index] * 1000);
	//	}
	//}
	//outds1->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, xsize, ysize, outb1data, xsize, ysize, GDT_Int16, 0, 0, 0);

GDALDataset* dst3 = (GDALDataset*)GDALOpen(b3Path.c_str(), GA_ReadOnly);
	GDALDataset* dst3Calibri = (GDALDataset*)GDALOpen(b3CalibriPath.c_str(), GA_ReadOnly);
	//dst3->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, originXsize, originYsize, origionData, originXsize, originYsize, GDT_UInt16, 0, 0);
	//dst3Calibri->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, calibriDataXsize, calibriDataYsize, calibridata, calibriDataXsize, calibriDataYsize, GDT_Float32, 0, 0);
	//for (int i = 0; i < ysize; i++)
	//{
	//	for (int j = 0; j < xsize; j++)
	//	{
	//		int line = lutline[i * xsize + j];
	//		int col = lutcol[i * xsize + j];
	//		//LatLonToRC(LonMin + step * j, LatMax - i * step, col, line);
	//		//outb1data[i * xsize + j] = origionData[originXsize * line + col];
	//		int index = origionData[originXsize * line + col];
	//		outb1data[i * xsize + j] = short(calibridata[index] * 1000);
	//	}
	//}
	//outds1->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, xsize, ysize, outb1data, xsize, ysize, GDT_Int16, 0, 0, 0);

//outds1->FlushCache();

//GDALClose(outds1);
	GDALDriver* driverMEM = GetGDALDriverManager()->GetDriverByName("MEM");
	GDALDriver* driverGTIFF = GetGDALDriverManager()->GetDriverByName("GTIFF");

double nReslution = 1000;
	double beginLineNum = 0;
	double geoTransform[6];
	geoTransform[0] = -5496000;
	geoTransform[1] = nReslution;
	geoTransform[2] = 0;
	geoTransform[3] = 5496000 - beginLineNum * nReslution;
	geoTransform[4] = 0;
	geoTransform[5] = -nReslution;
	dst1->SetGeoTransform(geoTransform);
	dst2->SetGeoTransform(geoTransform);
	dst3->SetGeoTransform(geoTransform);

GDALDataset* mergeDst = driverGTIFF->Create("F:\\testimgs\\rasterfiles\\FY4BProj.tif", 16236, 16205, 3, GDT_UInt16, NULL);

OGRSpatialReference* spatial = new OGRSpatialReference();
	spatial->importFromProj4("+proj=geos +no_defs +a=6378137.0 +b=6356752.3 +lon_0=133.0 +h=35785863");
	char* wkt = NULL;
	spatial->exportToWkt(&wkt);
	//std::cout << wkt << std::endl;

char** papszOptions = NULL;
	papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", wkt);
	papszOptions = CSLSetNameValue(papszOptions, "DST_SRS", wktName);
	//papszOptions = CSLSetNameValue(papszOptions, "SRC_METHOD", "NO_GEOTRANSFORM");

//void* hTransformArg = GDALCreateGenImgProjTransformer2(dst1, NULL, papszOptions);
	//void* hTransformArg = GDALCreateGenImgProjTransformer(dst3, wkt, NULL, wktName, FALSE, 0, 1);;
	double adfDstGeoTransform[6];
	//int nPixels = 0, nLines = 0;
	//GDALSuggestedWarpOutput(dst3, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines);
	adfDstGeoTransform[0] = 51.91;
	adfDstGeoTransform[1] = 0.01;
	adfDstGeoTransform[2] = 0;
	adfDstGeoTransform[3] = 81.05;
	adfDstGeoTransform[4] = 0;
	adfDstGeoTransform[5] = -0.01;
	//GDALDataset* outDatasetRGBproj = driverGTIFF->Create("F:\\testimgs\\rasterfiles\\FY4BProj.tif", nPixels, nLines, 1, GDT_UInt16, NULL);
	GDALDataset* outDataseMEM = driverMEM->Create("", 16236, 16205, 1, GDT_UInt16, NULL);
	outDataseMEM->SetProjection(wktName);
	outDataseMEM->SetGeoTransform(adfDstGeoTransform);
	warpdata(dst1, outDataseMEM);
	unsigned short* dataout = new unsigned short[16236 * 16205];
	outDataseMEM->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, 16236, 16205, dataout, 16236, 16205, GDT_UInt16, 0, 0);
	mergeDst->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, 16236, 16205, dataout, 16236, 16205, GDT_UInt16, 0, 0);
	warpdata(dst2, outDataseMEM);
	outDataseMEM->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, 16236, 16205, dataout, 16236, 16205, GDT_UInt16, 0, 0);
	mergeDst->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, 16236, 16205, dataout, 16236, 16205, GDT_UInt16, 0, 0);
	warpdata(dst3, outDataseMEM);
	outDataseMEM->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, 16236, 16205, dataout, 16236, 16205, GDT_UInt16, 0, 0);
	mergeDst->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, 16236, 16205, dataout, 16236, 16205, GDT_UInt16, 0, 0);
	//driverGTIFF->CreateCopy(, mergeDst, FALSE, 0, GDALTermProgress, 0);

mergeDst->SetGeoTransform(adfDstGeoTransform);
	mergeDst->SetProjection(wktName);

return 0;
}

```

扫描二维码推送至手机访问。

本文链接：https://lovedm.club/?id=154

分享给朋友：

返回列表

上一篇：C++ map使用

下一篇：linux g++编译gdal程序

个人学习笔记

FY4B warp

“FY4B warp” 的相关文章

C++ vector 初始化

C++ 计算点到空间直线的距离

C++ 函数指针

C++常量

GDAL构建Delaunay三角网

C++ 拷贝构造函数

鲁ICP备18000987号-1 鲁公网安备 37148202000241号

Powered By Z-BlogPHP. Theme by TOYEAN.