gdal 使用 RasterIO 进行栅格数据切片的时候,按照一副一副进行重采样效率太低,当等级比较高的时候,RasterIO 单瓦片速率都在17秒,这样我们就不能按照一张一张瓦片进行采样,而是降低RasterIO的使用次数,按照行列将数据全部读取到内存中,再进行数据输出就行
double TDRasterSplite::scale2Level(double scale) // //比例转换为层级
{
//实际距离
double mapdis = 0.0254 * ((double)tilesize / 96.0)*scale;
double x, y, x1, y1;
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, 0, 0, x, y);
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, tilesize, tilesize, x1, y1);
double srcscale = (x1 - x) / (0.0254 * ((double)tilesize / 96.0));
if (scale >= srcscale)
{
double level = mapdis / (x1 - x);
return fabs(level)*m_adfGeoTransform[1];
}
else
{
double s = mapdis / (x1 - x);
//x1 = x1 - mapdis; y1 = y1- mapdis;
//int c, r;
//this->LatLon2ColRow(adfGeoTransform, x1, y1, c, r);
//double level = (tilesize - r) / (double)tilesize;// -adfGeoTransform[1];
return fabs(s)*m_adfGeoTransform[1];
}
return 1;
}
bool TDRasterSplite::extractSpliteInfoByScale(std::string srcimagename, int tilesize, double scale, int& maxX, int& maxY)
{
GDALAllRegister();
m_srcimagename = srcimagename;
getOGRGeoTransform(); // 获取图层参数信息
this->tilesize = tilesize;
double level = scale2Level(scale);
maxX = 1.0*m_srcWidth / tilesize / level + 1;
maxY = 1.0*m_srcHeight / tilesize / level + 1;
return true;
}
bool TDRasterSplite::processByscale(std::string srcimagename, std::string outpath, int tilesize, std::string driverName, double scale, int startIndex, int spliteSize, int threadSize)
{
//初始化成员变量
GDALAllRegister();
m_format = driverName;
m_srcimagename = srcimagename;
m_outputstr = outpath;
getOGRGeoTransform(); // 获取图层参数信息
this->tilesize = tilesize;
double level = scale2Level(scale);
int xNum = 1.0*m_srcWidth / tilesize / level + 1;
int yNum = 1.0*m_srcHeight / tilesize / level + 1;
int beginX = startIndex % xNum;
int beginY = startIndex / xNum;
int curSpliteSize = 0;
TDThreadPool::ThreadPool threadPool(1, threadSize);
for (int y = beginY; y < yNum; y++) //行数
{
for (int x = beginX; x < xNum; x++) //列数
{
threadPool.CommitOneTask([this, xNum, yNum, x, y, level]() {
spliteoperator(xNum, yNum, x, y, level, m_outputstr, m_format, m_srcimagename, m_bandcount);
});
curSpliteSize++;
if (curSpliteSize >= spliteSize) {
goto SPLITE_END;
}
}
}
SPLITE_END:
threadPool.StopThreadPool();
return true;
}
bool TDRasterSplite::processByscale(std::string srcimagename, std::string outpath, int tilesize, std::string formattype, double scale, int le)
{
//初始化成员变量
GDALAllRegister();
m_format = formattype;
m_srcimagename = srcimagename;
m_outputstr = outpath;
getOGRGeoTransform(); // 获取图层参数信息
this->tilesize = tilesize;
double level = scale2Level(scale);
int xNum = 1.0*m_srcWidth / tilesize / level + 1;
int yNum = 1.0*m_srcHeight / tilesize / level + 1;
//TDThreadPool::ThreadPool threadPool(1, 1);
//for (int y = 0; y < yNum; y++) //行数
//{
// for (int x = 0; x < xNum; x++) //列数
// {
// threadPool.CommitOneTask([this, xNum, yNum, x, y, level]() {
// spliteoperator(xNum,yNum,x,y,level,m_outputstr,m_format,m_srcimagename,m_bandcount);
// });
// }
//}
//threadPool.StopThreadPool();
//TDThreadPool::ThreadPool threadPool(1, 5);
//for (int x = 0; x < xNum; x++) //列数
//{
// threadPool.CommitOneTask([this, xNum, yNum, x, level]() {
// spliteCol(xNum, yNum, x, level, m_outputstr, m_format, m_srcimagename, m_bandcount);
// });
//}
//threadPool.StopThreadPool();
TDThreadPool::ThreadPool threadPool(1, 1);
int x = 0;
int count = 20;
while (x < xNum) //列数
{
int xCount = count;
if (x + xCount > xNum) {
xCount = xNum - x;
}
threadPool.CommitOneTask([this, xNum, yNum, x, xCount, level]() {
spliteColNum(xNum, yNum, x, xCount,level, m_outputstr, m_format, m_srcimagename, m_bandcount);
});
x += xCount;
}
threadPool.StopThreadPool();
return true;
}
std::string TDRasterSplite::DecIntToHexStr(long long num)
{
std::string str;
long long Temp = num / 16;
int left = num % 16;
if (Temp > 0)
str += DecIntToHexStr(Temp);
if (left < 10)
str += (left + '0');
else
str += ('A' + left - 10);
if (str.length() < 4)
{
if (str.length() == 3) str = "0" + str;
if (str.length() == 2) str = "00" + str;
if (str.length() == 1) str = "000" + str;
}
return str;
}
std::string TDRasterSplite::DecStrToHexStr(std::string str)
{
long long Dec = 0;
for (int i = 0; i < str.size(); ++i)
Dec = Dec * 10 + str[i] - '0';
return DecIntToHexStr(Dec);
}
double* TDRasterSplite::UpdateGeoTransfrom(double adfGeoTransform[], double minx, double maxy, double l)
{
double resadfGeoTransform[6];
resadfGeoTransform[0] = minx;
resadfGeoTransform[1] = l * adfGeoTransform[1];
resadfGeoTransform[2] = 0;
resadfGeoTransform[3] = maxy;
resadfGeoTransform[4] = 0;
resadfGeoTransform[5] = l * adfGeoTransform[5];
return resadfGeoTransform;
}
bool TDRasterSplite::create_copy(GDALDataset* dataset, std::string tiffile, GDALDriver* jpg_driver)
{
if (!dataset)
return false;
std::string jpgfile = tiffile;
if (!jpg_driver)
return false;
GDALDataset*jpg_dataset = nullptr;
jpg_dataset = jpg_driver->CreateCopy(jpgfile.c_str(), dataset, true, 0, 0, 0);
if (!jpg_dataset)
return false;
GDALClose(jpg_dataset);
return true;
}
void TDRasterSplite::spliteColNum(int xNum, int yNum, int x, int xCount, double level, std::string filepath, std::string formattype, std::string srcimagename, int bandCount)
{
std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
auto srcDataset = (GDALDataset*)GDALOpen(srcimagename.c_str(), GA_ReadOnly); //获取原始数据句柄
if (nullptr == srcDataset)
{
return;
}
auto armDriver = GetGDALDriverManager()->GetDriverByName(formattype.c_str());
//GDALRasterIOExtraArg exterArg;
//INIT_RASTERIO_EXTRA_ARG(exterArg);
//exterArg.eResampleAlg = GDALRIOResampleAlg::GRIORA_Bilinear;
//读取影响的位置
double startX = x * tilesize *level;
double startY = 0;
double sw = (m_srcWidth - (xNum - 1)*tilesize * level) / (tilesize * level);
double w = m_srcWidth - (xNum - 1)*tilesize * level;
double sh = (m_srcHeight - (yNum - 1)*tilesize * level) / (tilesize * level);
double h = m_srcHeight - (yNum - 1)*tilesize * level;
//
int storeWidth = 0;
//读取影像的宽度
int readWidth = 0;
if (x + xCount == xNum) {
readWidth = tilesize*(xCount - 1)*level + w;
storeWidth = tilesize * (xCount - 1) + tilesize * sw;
}
else {
readWidth = tilesize*level* xCount;
storeWidth = tilesize * xCount;
}
//读取影像的高度
int readHeight = m_srcHeight;
int storeHeight = tilesize * (yNum - 1) + tilesize * sh;
//开辟空间
size_t imgBufNum = storeWidth * storeHeight * bandCount;
GByte *imgBuf = new GByte[imgBufNum];
memset(imgBuf,0, imgBufNum);
//读取N列瓦片影像
srcDataset->RasterIO(GF_Read, startX, startY, readWidth, readHeight, imgBuf, storeWidth, storeHeight,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
std::cout << "TotalTime : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;
TDThreadPool::ThreadPool threadPool(1, 20);
for (int y = 0; y < yNum; ++y) {
int tileSizeX = tilesize;
int tileSizeY = tilesize;
if ((yNum - 1) == y) {
tileSizeY = tilesize * sh;
}
for (int i = 0; i < xCount; ++i) {
int curX = x + i;
if ((xNum - 1) == curX) {
tileSizeX = tilesize * sw;
}
std::string target_name = filepath + "/R" + DecIntToHexStr(y) + "C" + DecIntToHexStr(curX) + TDAPICommon::GetSuffixName(srcimagename); // 配置输出参数
threadPool.CommitOneTask([this, curX,i,y, level,tileSizeX, tileSizeY, storeWidth, storeHeight, armDriver,target_name, bandCount, srcDataset, imgBuf]()
{
GDALDataset* poDstDS;
poDstDS = armDriver->Create(target_name.c_str(), tilesize, tilesize, bandCount, GDT_Byte, nullptr);
//设置放射变换参数
CPLErr error = poDstDS->SetProjection(srcDataset->GetProjectionRef());
double startX = curX * tilesize*level;
double startY = y * tilesize*level;
double xx, yy;
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, startX, startY, xx, yy);
poDstDS->SetGeoTransform(UpdateGeoTransfrom(m_adfGeoTransform, xx, yy, level)); //设置图像参数
GByte *pStart = new GByte[tileSizeX*tileSizeY*bandCount];
memset(pStart, 0, tileSizeX*tileSizeY*bandCount);
//计算原始影像一行瓦片的长度
int rowTileWidth = storeWidth * tilesize;
//计算N行瓦片的长度
int nRowTileWidth = rowTileWidth * y;
for (int iBand = 0;iBand < bandCount;++iBand) {
//计算数据瓦片的波段起始位置
int dstStart = tileSizeX * tileSizeY * iBand;
//计算读取的原始影像的影像的起始位置
int srcStart = storeWidth * storeHeight * iBand;
for (int j = 0; j < tileSizeY; ++j) {
memcpy(pStart + dstStart + j * tileSizeX, imgBuf + srcStart + nRowTileWidth + j* storeWidth + i * tilesize, tileSizeX);
}
}
poDstDS->RasterIO(GF_Write, 0, 0, tileSizeX, tileSizeY, pStart, tileSizeX, tileSizeY,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
GDALClose(poDstDS);
poDstDS = NULL;
delete[] pStart;
});
}
}
threadPool.StopThreadPool();
delete[] imgBuf;
imgBuf = nullptr;
GDALClose(srcDataset);
srcDataset = NULL;
/*std::chrono::steady_clock::time_point */end = std::chrono::steady_clock::now();
std::cout << "TotalTime : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;
}
void TDRasterSplite::spliteCol(int xNum, int yNum, int x, double level, std::string filepath, std::string formattype, std::string srcimagename, int bandCount)
{
std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
auto srcDataset = (GDALDataset*)GDALOpen(srcimagename.c_str(), GA_ReadOnly); //获取原始数据句柄
if (nullptr == srcDataset)
{
return;
}
auto armDriver = GetGDALDriverManager()->GetDriverByName(formattype.c_str());
GDALRasterIOExtraArg exterArg;
INIT_RASTERIO_EXTRA_ARG(exterArg);
exterArg.eResampleAlg = GDALRIOResampleAlg::GRIORA_NearestNeighbour;
//读取影响的位置
int startX = x * tilesize *level;
int startY = 0;
double sw = (m_srcWidth - (xNum - 1)*tilesize * level) / (double)tilesize / (double)level;
int w = m_srcWidth - (xNum - 1)*tilesize * level;
double sh = (m_srcHeight - (yNum - 1)*tilesize * level) / (double)tilesize / (double)level;
int h = m_srcHeight - (yNum - 1)*tilesize * level;
//存储影像的宽度
int storeWidth = tilesize;
//读取影像的宽度
int readWidth = tilesize * level;
if ((xNum - 1) == x) {
readWidth = w;
storeWidth = tilesize*sw;
}
//读取影像的高度
size_t readHeight = tilesize * (yNum - 1)*level + h;
size_t stortHeight = tilesize * (yNum - 1) + tilesize * sh;
//开辟空间
size_t imgBufNum = (size_t)(storeWidth*stortHeight) * bandCount;
GByte *imgBuf = new GByte[imgBufNum];
//读取一列瓦片影像
srcDataset->RasterIO(GF_Read, startX, startY, readWidth, readHeight, imgBuf, storeWidth, stortHeight,
GDT_Byte, bandCount, nullptr, 0, 0, 0, &exterArg);
for (int y = 0; y < yNum; ++y) {
int tileSizeX = tilesize;
int tileSizeY = tilesize;
if ((yNum - 1) == y) {
tileSizeY = tilesize * sh;
}
if ((xNum - 1) == x) {
tileSizeX = storeWidth;
}
std::string target_name1 = filepath + "/R" + DecIntToHexStr(y) + "C" + DecIntToHexStr(x) + TDAPICommon::GetSuffixName(srcimagename); // 配置输出参数
GDALDataset* poDstDS;
poDstDS = armDriver->Create(target_name1.c_str(), tilesize, tilesize, bandCount, GDT_Byte, nullptr);
//设置放射变换参数
CPLErr error = poDstDS->SetProjection(srcDataset->GetProjectionRef());
startY = y * tilesize *level;
double xx, yy, xx1, yy1;
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, startX, startY, xx, yy);
poDstDS->SetGeoTransform(UpdateGeoTransfrom(m_adfGeoTransform, xx, yy, level)); //设置图像参数
GByte *pStart = imgBuf + (tileSizeX*tileSizeY)*sizeof(GByte);
poDstDS->RasterIO(GF_Write, 0, 0, tileSizeX, tileSizeY, pStart, tileSizeX, tileSizeY,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
GDALClose(poDstDS);
poDstDS = NULL;
}
delete[] imgBuf;
imgBuf = nullptr;
GDALClose(srcDataset);
srcDataset = NULL;
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
std::cout << "TotalTime : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;
}
void TDRasterSplite::spliteoperator(int xNum, int yNum, int x, int y, double level, std::string filepath,std::string formattype, std::string srcimagename,int bandCount)
{
//std::string target_name = filepath + "/R" + DecIntToHexStr(y) + "C" + DecIntToHexStr(x) + ".tif";
//std::cout << target_name << "x:"<<x<< "y:"<<y<<std::endl;
std::string target_name1 = filepath + "/R" + DecIntToHexStr(y) + "C" + DecIntToHexStr(x) + TDAPICommon::GetSuffixName(srcimagename); // 配置输出参数
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
auto srcDataset = (GDALDataset*)GDALOpen(srcimagename.c_str(), GA_ReadOnly); //获取原始数据句柄
if (nullptr == srcDataset)
{
return;
}
auto armDriver = GetGDALDriverManager()->GetDriverByName(formattype.c_str());
auto poDriver = GetGDALDriverManager()->GetDriverByName("MEM");
int startX = x * tilesize *level;
int startY = y * tilesize*level;
GDALDataset* poDstDS;
poDstDS = armDriver->Create(target_name1.c_str(), tilesize, tilesize, bandCount, GDT_Byte, nullptr);
//设置放射变换参数
CPLErr error = poDstDS->SetProjection(srcDataset->GetProjectionRef());
double xx, yy, xx1, yy1;
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, startX, startY, xx, yy);
poDstDS->SetGeoTransform(UpdateGeoTransfrom(m_adfGeoTransform, xx, yy, level)); //设置图像参数
//GDALRasterIOExtraArg exterArg;
//INIT_RASTERIO_EXTRA_ARG(exterArg);
//exterArg.eResampleAlg = GDALRIOResampleAlg::GRIORA_Bilinear;
if (y == yNum - 1 || x == xNum - 1) // 进行重采样
{
double sw = 1; double sh = 1;
int w = tilesize * level;
int h = tilesize * level;
if (x == xNum - 1)
{
sw = (m_srcWidth - (xNum - 1)*tilesize * level) / (double)tilesize / (double)level;
w = m_srcWidth - (xNum - 1)*tilesize * level;
}
if (y == yNum - 1)
{
sh = (m_srcHeight - (yNum - 1)*tilesize * level) / (double)tilesize / (double)level;
h = m_srcHeight - (yNum - 1)*tilesize * level;
}
size_t imgBufNum = (size_t)tilesize*sw * tilesize*sh * bandCount;
GByte *imgBuf = new GByte[imgBufNum];
srcDataset->RasterIO(GF_Read, startX, startY, w, h, imgBuf, tilesize*sw, tilesize*sh,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
poDstDS->RasterIO(GF_Write, 0, 0, tilesize*sw, tilesize*sh, imgBuf, tilesize*sw, tilesize*sh,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
//create_copy(poDstDS, target_name1, armDriver);
//释放资源
delete[] imgBuf;
imgBuf = nullptr;
}
else
{
size_t imgBufNum = (size_t)tilesize * tilesize * bandCount;
GByte *imgBuf = new GByte[imgBufNum];
std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
srcDataset->RasterIO(GF_Read, startX, startY, tilesize*level, tilesize*level, imgBuf, tilesize, tilesize,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
std::cout << "TotalTime : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s";
poDstDS->RasterIO(GF_Write, 0, 0, tilesize, tilesize, imgBuf, tilesize, tilesize,
GDT_Byte, bandCount, nullptr, 0, 0, 0);
//create_copy(poDstDS, target_name1, armDriver);
//释放资源
delete[] imgBuf;
imgBuf = nullptr;
}
GDALClose(poDstDS);
poDstDS = NULL;
GDALClose(srcDataset);
srcDataset = NULL;
}
void TDRasterSplite::getOGRGeoTransform()
{
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
auto srcDataset = (GDALDataset*)GDALOpen(m_srcimagename.c_str(), GA_ReadOnly);
if (srcDataset == nullptr)
return;
srcDataset->GetGeoTransform(m_adfGeoTransform);
m_srcWidth = srcDataset->GetRasterXSize();
m_srcHeight = srcDataset->GetRasterYSize();
m_bandcount = srcDataset->GetRasterCount();
GDALDriver* poDriver = srcDataset->GetDriver();
char** papszMetadata;
papszMetadata = poDriver->GetMetadata();
m_drivateflags = CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE);
GDALClose(srcDataset);
srcDataset = NULL;
}
double TDRasterSplite::level2scale(double level) //计算当前层级比例尺
{
//return 1.0 / (0.0254 / (2 * 96));
//地图距离
double mapdis = 0.0254 * ((double)tilesize / 96.0);
//实际距离
double x, y, x1, y1;
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, 0, 0, x, y);
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, tilesize*level, tilesize*level, x1, y1);
return (x1 - x) / mapdis;
}
bool TDRasterSplite::process(LayerPathInfoV1 originalPathInfo, double grid_level, std::string pszFormat,std::string outpath,std::map<std::string, std::tuple<double, double, double, double>> splitinformation, std::map<std::string, std::string> &objpathinfo)
{
// 获取配置参数
TDRasterCommon obj;
obj.process(originalPathInfo.filePath);
/*GT(0)和GT(3)是第一组,表示图像左上角的地理坐标;
GT(1)和GT(5)是第二组,表示图像横向和纵向的分辨率(一般这两者的值相等,符号相反,横向分辨率为正数,纵向分辨率为负数);
GT(2)和GT(4)是第三组,表示图像旋转系数,对于一般图像来说,这两个值都为0*/
m_adfGeoTransform[0] = obj.m_adfGeoTransform[0];
m_adfGeoTransform[1] = obj.m_adfGeoTransform[1];
m_adfGeoTransform[2] = obj.m_adfGeoTransform[2];
m_adfGeoTransform[3] = obj.m_adfGeoTransform[3];
m_adfGeoTransform[4] = obj.m_adfGeoTransform[4];
m_adfGeoTransform[5] = obj.m_adfGeoTransform[5];
m_bandcount = obj.m_bandcount;
m_srcHeight = obj.m_srcHeight;
m_srcWidth = obj.m_srcWidth;
OGREnvelope allogrenvelope;
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, 0, 0, allogrenvelope.MinX,allogrenvelope.MinY);
TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, m_srcHeight,m_srcWidth,allogrenvelope.MaxX, allogrenvelope.MaxY);
//auto ogrsrs = obj.m_ogr_srs;
//TDRasterCommon::CoordinateTransformation(ogrsrs, allogrenvelope.MinX, allogrenvelope.MinY, allogrenvelope.MinX, allogrenvelope.MinY);
//TDRasterCommon::CoordinateTransformation(ogrsrs, allogrenvelope.MaxX, allogrenvelope.MaxY, allogrenvelope.MaxX, allogrenvelope.MaxY);
if (allogrenvelope.MinX > allogrenvelope.MaxX)
{
auto temp = allogrenvelope.MaxX;
allogrenvelope.MaxX = allogrenvelope.MinX;
allogrenvelope.MinX = temp;
}
if (allogrenvelope.MinY > allogrenvelope.MaxY)
{
auto temp = allogrenvelope.MaxY;
allogrenvelope.MaxY = allogrenvelope.MinY;
allogrenvelope.MinY = temp;
}
std::for_each(splitinformation.begin(), splitinformation.end(), [&](std::map<std::string, std::tuple<double, double, double, double>>::value_type &it) {
OGREnvelope tempnvelope, ogrenvelope;
ogrenvelope = allogrenvelope;
tempnvelope.MinX = std::get<0>(it.second);
tempnvelope.MinY = std::get<1>(it.second);
tempnvelope.MaxX = std::get<2>(it.second);
tempnvelope.MaxY = std::get<3>(it.second);
if (ogrenvelope.Intersects(tempnvelope) || ogrenvelope.Contains(tempnvelope))
{
ogrenvelope.Intersect(tempnvelope);
m_mapOGREnvelope[it.first] = ogrenvelope;
m_pathinfo[it.first];
}
if (tempnvelope.Intersects(ogrenvelope) || tempnvelope.Contains(ogrenvelope))
{
ogrenvelope.Intersect(tempnvelope);
m_mapOGREnvelope[it.first] = ogrenvelope;
m_pathinfo[it.first];
}
});
TDThreadPool::ThreadPool threadPool(10, 10);
std::for_each(m_mapOGREnvelope.begin(), m_mapOGREnvelope.end(), [&](std::map<std::string, OGREnvelope>::value_type it) {
auto filename = TDAPICommon::GetUniqueDirName();
std::string srcFile;
if (0 == pszFormat.compare("img"))
{
srcFile = TDAPICommon::JoinPath(outpath, filename + ".img");
}
if (0 == pszFormat.compare("GTiff"))
{
srcFile = TDAPICommon::JoinPath(outpath, filename + ".tif");
}
double fResX = grid_level / ABS(m_adfGeoTransform[1]);
double fResY = grid_level / ABS(m_adfGeoTransform[5]);
auto filePath = originalPathInfo.filePath;
OGREnvelope params = it.second;
auto pszFormatParams = pszFormat;
//auto temp = new OGRSpatialReference();
//temp->importFromWkt(R"~(PROJCS["CGCS2000_3_Degree_GK_Zone_45",GEOGCS["GCS_China_Geodetic_Coordinate_System_2000",DATUM["D_China_2000",SPHEROID["CGCS2000",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",45500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",135.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG",4533]])~");
//TDRasterCommon::CoordinateTransformation(temp, params.MinX, params.MinY, params.MinX, params.MinY);
//TDRasterCommon::CoordinateTransformation(temp, params.MaxX, params.MaxY, params.MaxX, params.MaxY);
threadPool.CommitOneTask([this,filePath,srcFile, fResX, fResY, params, pszFormatParams]() {
CutImageByAOIGDAL(filePath, srcFile, fResX, fResY, GDALResampleAlg(0), { params.MinX,params.MinY,params.MaxX, params.MaxY}, NULL, 0, pszFormatParams);
});
m_pathinfo[it.first] = srcFile;
});
threadPool.StopThreadPool();
objpathinfo = m_pathinfo;
//获取需要切分的四至点信息
return true;
}
/** 重采样函数(GDAL)
* @param pszSrcFile 输入文件的路径
* @param pszOutFile 写入的结果图像的路径
* @param fResX X转换采样比,默认大小为1.0,大于1图像变大,小于1表示图像缩小
* @param fResY Y转换采样比,默认大小为1.0
* @param nResampleMode 采样模式,有五种,具体参见GDALResampleAlg定义,默认为双线性内插
* @param pExtent 采样范围,为NULL表示计算全图
* @param pBandIndex 指定的采样波段序号,为NULL表示采样全部波段
* @param pBandCount 采样的波段个数,同pBandIndex一同使用,表示采样波段的个数
* @param pszFormat 写入的结果图像的格式
* @return 成功返回0,否则为其他值
*/
bool TDRasterSplite::CutImageByAOIGDAL(std::string pszSrcFile, std::string pszOutFile, float fResX, float fResY, GDALResampleAlg nResampleMode,
std::vector<double> vecdouble, int *pBandIndex, int pBandCount, std::string pszFormat)
{
if (vecdouble.size() < 4)
{
return false;
}
OGREnvelope pEnvelope;
pEnvelope.MinX = vecdouble[0];
pEnvelope.MinY = vecdouble[1];
pEnvelope.MaxX = vecdouble[2];
pEnvelope.MaxY = vecdouble[3];
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile.c_str(), GA_ReadOnly);
if (pDSrc == NULL)
{
return false;
}
GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat.c_str());
if (pDriver == NULL)
{
GDALClose((GDALDatasetH)pDSrc);
return false;
}
int iBandCount = pDSrc->GetRasterCount();
std::string strWkt = pDSrc->GetProjectionRef();
GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
double dGeoTrans[6] = { 0 };
double adfGeoTransform[6] = { 0 };
pDSrc->GetGeoTransform(dGeoTrans);
pDSrc->GetGeoTransform(adfGeoTransform);
int iNewBandCount = iBandCount;
if (pBandIndex != NULL && pBandCount != NULL)
{
int iMaxBandIndex = pBandIndex[0];
for (int i = 1; i < pBandCount; i++)
{
if (iMaxBandIndex < pBandIndex[i])
iMaxBandIndex = pBandIndex[i];
}
if (iMaxBandIndex > iBandCount)
{
GDALClose((GDALDatasetH)pDSrc);
return false;
}
iNewBandCount = pBandCount;
}
double dPrj[4] = { 0 };
int iNewWidth = 0, iNewHeigh = 0;
int iBeginRow = 0, iBeginCol = 0;
dGeoTrans[0] = pEnvelope.MinX;
dGeoTrans[3] = pEnvelope.MaxY;
dGeoTrans[1] = dGeoTrans[1] * fResX;
dGeoTrans[5] = dGeoTrans[5] * fResY;
iNewWidth = static_cast<int>((pEnvelope.MaxX - pEnvelope.MinX) / ABS(dGeoTrans[1]) + 0.5);
iNewHeigh = static_cast<int>((pEnvelope.MaxY - pEnvelope.MinY) / ABS(dGeoTrans[5]) + 0.5);
TDRasterCommon::Projection2ImageRowCol(adfGeoTransform, dGeoTrans[0], dGeoTrans[3], iBeginCol, iBeginRow);
TDRasterCommon::ImageRowCol2Projection(adfGeoTransform, iBeginCol, iBeginRow, dGeoTrans[0], dGeoTrans[3]);
GDALDataset *pDDst = pDriver->Create(pszOutFile.c_str(), iNewWidth, iNewHeigh, iNewBandCount, dataType, NULL);
if (pDDst == NULL)
{
//GDALClose((GDALDatasetH)pDSrc);
return false;
}
pDDst->SetProjection(strWkt.c_str());
pDDst->SetGeoTransform(dGeoTrans);
GDALResampleAlg eResample = (GDALResampleAlg)nResampleMode;
int *pSrcBand = NULL;
int *pDstBand = NULL;
int iBandSize = 0;
if (pBandIndex != NULL && pBandCount != NULL)
{
iBandSize = pBandCount;
pSrcBand = new int[iBandSize];
pDstBand = new int[iBandSize];
for (int i = 0; i < iBandSize; i++)
{
pSrcBand[i] = pBandIndex[i];
pDstBand[i] = i + 1;
}
}
else
{
iBandSize = iBandCount;
pSrcBand = new int[iBandSize];
pDstBand = new int[iBandSize];
for (int i = 0; i < iBandSize; i++)
{
pSrcBand[i] = i + 1;
pDstBand[i] = i + 1;
}
}
void *hTransformArg = NULL, *hGenImgPrjArg = NULL;
hTransformArg = hGenImgPrjArg = GDALCreateGenImgProjTransformer2((GDALDatasetH)pDSrc, (GDALDatasetH)pDDst, NULL);
if (hTransformArg == NULL)
{
GDALClose((GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)pDDst);
return false;
}
GDALTransformerFunc pFnTransformer = GDALGenImgProjTransform;
GDALWarpOptions *psWo = GDALCreateWarpOptions();
psWo->papszWarpOptions = CSLDuplicate(NULL);
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
psWo->pfnTransformer = pFnTransformer;
psWo->pTransformerArg = hTransformArg;
//psWo->pfnProgress = GDALTermProgress;
psWo->nBandCount = iNewBandCount;
psWo->panSrcBands = (int *)CPLMalloc(iNewBandCount * sizeof(int));
psWo->panDstBands = (int *)CPLMalloc(iNewBandCount * sizeof(int));
for (int i = 0; i < iNewBandCount; i++)
{
psWo->panSrcBands[i] = pSrcBand[i];
psWo->panDstBands[i] = pDstBand[i];
}
delete[] pSrcBand;
delete[] pDstBand;
GDALWarpOperation oWo;
if (oWo.Initialize(psWo) != CE_None)
{
GDALClose((GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)pDDst);
return false;
}
oWo.ChunkAndWarpImage(0, 0, iNewWidth, iNewHeigh);
GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)pDDst);
return true;
}