GDAL 根据图版范围对影像进行提取 (C++篇)


bool ImageExtract::CutImageByFeature(string shpfile, string field1, string field2, string savePath)
{std::string message;//判断影像if (imageDataset == nullptr) {return false;}savePath += "/";const char* image_projection = imageDataset->GetProjectionRef();OGRSpatialReference image_srs;image_srs.importFromWkt(image_projection);if ((image_projection == nullptr) || !image_srs.IsProjected()){//非投影坐标系统return false;}GDALDataset *poDS = nullptr;OGRLayer *poLayer = nullptr;LayerPathInfoV1 feature_info;feature_info.drivername = "ESRI Shapefile";feature_info.filename = filename(shpfile);feature_info.filePath = shpfile;open_layer(feature_info, poDS, poLayer);if (!poDS || !poLayer) {return false;}OGRSpatialReference* feature_srs = poLayer->GetSpatialRef();if ((feature_srs == nullptr) /*|| !feature_srs->IsProjected()*/){//非投影坐标系统return false;}if (!feature_srs->IsProjected()){QString newshp = QString::fromStdString(shpfile);QFileInfo fileinfo(newshp);QString newname = fileinfo.path() + "/" + fileinfo.fileName().split(".")[0] + "_new.shp";//地理坐标转平面坐标if (TransGeo2Project(poLayer, newname.toStdString(), &image_srs)) {is_shp_created_ = true;new_shpfile = newname.toUtf8().toStdString();//关闭好的if (poDS) {GDALClose(poDS);poDS = nullptr;poLayer = nullptr;feature_srs = nullptr;}feature_info.filename = filename(new_shpfile);feature_info.filePath = new_shpfile;open_layer(feature_info, poDS, poLayer);if (!poDS || !poLayer) {return false;}feature_srs = poLayer->GetSpatialRef();if ((feature_srs == nullptr)){return false;}}else {return false;}}//判断俩者坐标系相同if (!IsSameSrs(feature_srs, &image_srs)) {return false;}//判断影像与矢量相交double image_geomtransform[6] = { 0 };imageDataset->GetGeoTransform(image_geomtransform);int image_size_x = imageDataset->GetRasterXSize();int image_size_y = imageDataset->GetRasterYSize();OGREnvelope image_envelope;OGRGeometry* image_geometry = get_iamge_envlope(image_geomtransform,image_size_x, image_size_y, image_envelope);if (nullptr == image_geometry) {return false;}OGREnvelope *LayerEnvelope = new OGREnvelope();poLayer->GetExtent(LayerEnvelope);if (!LayerEnvelope->Intersects(image_envelope)) {return false;}char* wkt = nullptr;poLayer->GetSpatialRef()->exportToWkt(&wkt);//判别影像分辨率int Xpixels = imageDataset->GetRasterXSize();int Ypixels = imageDataset->GetRasterYSize();double xmin, ymin, xmax, ymax;ColRow2LatLon(adfGeoTransform, 0, 0, xmin, ymax);ColRow2LatLon(adfGeoTransform, Xpixels, Ypixels, xmax, ymin);OGRFeature *poFeature;poLayer->ResetReading();//int ccc = 0;int inext = 0, itotal = poLayer->GetFeatureCount();int iinternal = itotal / 20, ishow = 0, icur = 0;while ((poFeature = poLayer->GetNextFeature()) != NULL){icur++;if (icur > inext && ishow <= 100) {inext += iinternal;ishow += 5;}OGREnvelope FeatureEnvelop;if(poFeature->GetGeometryRef() == nullptr)continue;poFeature->GetGeometryRef()->getEnvelope(&FeatureEnvelop);///ccc++;double s1 = (FeatureEnvelop.MaxX - FeatureEnvelop.MinX) / image_geomtransform[1] / defaultwid;double s2 = (FeatureEnvelop.MaxY - FeatureEnvelop.MinY) / image_geomtransform[1] / defaultwid;double s = (s1 > s2) ? s1 : s2;//合并名称string ccc = poFeature->GetFieldAsString(field1.c_str());if (ccc.length() > 6)ccc = ccc.substr(0, 6);std::string tbbh_value = poFeature->GetFieldAsString(field2.c_str());ccc = ccc + "_" + tbbh_value;if (s < 1){//计算地理坐标double w = defaultwid * (xmax - xmin) / Xpixels - FeatureEnvelop.MaxX + FeatureEnvelop.MinX;double h = defaultheight * (ymax - ymin) / Ypixels - FeatureEnvelop.MaxY + FeatureEnvelop.MinY;double lx = FeatureEnvelop.MinX - w / 2;double rx = FeatureEnvelop.MaxX + w / 2;double ty = FeatureEnvelop.MinY - h / 2;double by = FeatureEnvelop.MaxY + h / 2;int c0, r0, c1, r1;LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);int blockex_width = fabs(r1 - r0); int  blockex_height = fabs(c1 - c0);if (blockex_width <= 0 || blockex_height <= 0)continue;double geox = min(lx, rx), geoy = max(by, ty);//处理边界问题if ((r0 < 0 || c0 < 0)&& image_envelope.Contains(FeatureEnvelop)) {double w = FeatureEnvelop.MaxX - FeatureEnvelop.MinX;double h = FeatureEnvelop.MaxY - FeatureEnvelop.MinY;double extent = 0;if (w > h)extent = w * extend_factor;else extent = h * extend_factor;lx = FeatureEnvelop.MinX - extent;rx = FeatureEnvelop.MaxX + extent;ty = FeatureEnvelop.MinY - extent;by = FeatureEnvelop.MaxY + extent;              LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);//再次保证切片影像if (r0 < 0 || c0 < 0) {extent = w * 0;lx = FeatureEnvelop.MinX - extent;rx = FeatureEnvelop.MaxX + extent;ty = FeatureEnvelop.MinY - extent;by = FeatureEnvelop.MaxY + extent;LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);}blockex_width = fabs(r1 - r0); blockex_height = fabs(c1 - c0);//宽度长度为0返回if (blockex_width <= 0 || blockex_height <= 0)continue;geox = min(lx, rx), geoy = max(by, ty);}if (blockex_width > (Xpixels - r0))blockex_width = Xpixels - r0;if (blockex_height > (Ypixels - c0))blockex_height = Ypixels - c0;CreateSubimg(poFeature->GetGeometryRef(), r0, c0, blockex_width, blockex_height, geox, geoy, wkt, ccc, savePath);}else if (s >= 1){//计算金字塔级别double w = FeatureEnvelop.MaxX - FeatureEnvelop.MinX;double h = FeatureEnvelop.MaxY - FeatureEnvelop.MinY;double extent = 0;if (w > h){extent = w * extend_factor;//左右下上double lx = (FeatureEnvelop.MaxX - FeatureEnvelop.MinX) / 2.0 + FeatureEnvelop.MinX - w / 2 - extent;double rx = (FeatureEnvelop.MaxX - FeatureEnvelop.MinX) / 2.0 + FeatureEnvelop.MinX + w / 2 + extent;double ty = (FeatureEnvelop.MaxX - FeatureEnvelop.MinX) / 2.0 + FeatureEnvelop.MinY - w / 2 - extent;double by = (FeatureEnvelop.MaxX - FeatureEnvelop.MinX) / 2.0 + FeatureEnvelop.MinY + w / 2 + extent;//c0:上行 r0:左列 c1:下行 r1:右列int c0, r0, c1, r1;LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);int blockex_width = fabs(r1 - r0); int  blockex_height = fabs(c1 - c0);if (blockex_width <= 0 || blockex_height <= 0)continue;//判断扩展大于影像情况并且图版包含或部分包含在影像内情况if ((r0 < 0 || c0 < 0)&& image_envelope.Contains(FeatureEnvelop)) {extent = w * 0;lx = FeatureEnvelop.MinX - extent;rx = FeatureEnvelop.MaxX + extent;ty = FeatureEnvelop.MinY - extent;by = FeatureEnvelop.MaxY + extent;//c0:上行 r0:左列 c1:下行 r1:右列LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);blockex_width = fabs(r1 - r0); blockex_height = fabs(c1 - c0);if (blockex_width <= 0 || blockex_height <= 0)continue;}else {//影像完全包含图版或者影像没有包含图版}if (blockex_width > (Xpixels - r0))blockex_width = Xpixels - r0;if (blockex_height > (Ypixels - c0))blockex_height = Ypixels - c0;CreateSubimg(poFeature->GetGeometryRef(), r0, c0, blockex_width, blockex_height, min(lx, rx), max(by, ty), wkt, ccc, savePath);}else{extent = h * extend_factor;//左右下上double lx = (FeatureEnvelop.MaxY - FeatureEnvelop.MinY) / 2.0 + FeatureEnvelop.MinX - h / 2 - extent;double rx = (FeatureEnvelop.MaxY - FeatureEnvelop.MinY) / 2.0 + FeatureEnvelop.MinX + h / 2 + extent;double ty = (FeatureEnvelop.MaxY - FeatureEnvelop.MinY) / 2.0 + FeatureEnvelop.MinY - h / 2 - extent;double by = (FeatureEnvelop.MaxY - FeatureEnvelop.MinY) / 2.0 + FeatureEnvelop.MinY + h / 2 + extent;//c0:上行 r0:左列 c1:下行 r1:右列int c0, r0, c1, r1;LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);int blockex_width = fabs(r1 - r0); int  blockex_height = fabs(c1 - c0);if(blockex_width <= 0 || blockex_height <= 0)continue;//判断扩展大于影像情况并且图版包含或部分包含在影像内情况if ((r0 < 0 || c0 < 0)&& image_envelope.Contains(FeatureEnvelop)) {extent = w * 0;lx = FeatureEnvelop.MinX - extent;rx = FeatureEnvelop.MaxX + extent;ty = FeatureEnvelop.MinY - extent;by = FeatureEnvelop.MaxY + extent;//c0:上行 r0:左列 c1:下行 r1:右列LatLon2ColRow(adfGeoTransform, lx, ty, c1, r0);LatLon2ColRow(adfGeoTransform, rx, by, c0, r1);blockex_width = fabs(r1 - r0);blockex_height = fabs(c1 - c0);if (blockex_width <= 0 || blockex_height <= 0)continue;}else {//影像完全包含图版或者影像没有包含图版}if (blockex_width > (Xpixels - r0))blockex_width = Xpixels - r0;if (blockex_height > (Ypixels - c0))blockex_height = Ypixels - c0;CreateSubimg(poFeature->GetGeometryRef(), r0, c0, blockex_width, blockex_height, min(lx,rx), max(by,ty), wkt, ccc, savePath);}}OGRFeature::DestroyFeature(poFeature);}if(poDS)GDALClose(poDS);return true;
}bool ImageExtract::CreateSubimg(OGRGeometry *geometry, int c0, int r0, int blockex_width, int blockex_height, double minx, double maxy, const char *wkt, string ccc, string SavePath)
{char* blockex_1 = nullptr, *blockex_2 = nullptr, *blockex_3 = nullptr;bool isblockex_1 = read_block(image_band1, c0, r0, blockex_width, blockex_height, blockex_1);bool isblockex_2 = read_block(image_band2, c0, r0, blockex_width, blockex_height, blockex_2);bool isblockex_3 = read_block(image_band3, c0, r0, blockex_width, blockex_height, blockex_3);GDALDataset* output_dataset = nullptr;string name = SavePath + ccc + ".tif";create_rasterband("GTiff", name, 3, blockex_width, blockex_height, output_dataset);resadfGeoTransform[0] = minx;resadfGeoTransform[1] = adfGeoTransform[1];resadfGeoTransform[2] = 0;resadfGeoTransform[3] = maxy;resadfGeoTransform[4] = 0;resadfGeoTransform[5] = adfGeoTransform[5];output_dataset->SetGeoTransform(resadfGeoTransform);if (output_dataset != NULL)output_dataset->SetProjection(wkt);GDALRasterBand* merge_band1 = output_dataset->GetRasterBand(1);GDALRasterBand* merge_band2 = output_dataset->GetRasterBand(2);GDALRasterBand* merge_band3 = output_dataset->GetRasterBand(3);if (!merge_band1 || !merge_band2 || !merge_band3) {//message = "Error : get tmp image band failure, feature id :" + icur;//ibreak = false;if (output_dataset)GDALClose(output_dataset);return false;}//DrawGeometry(geometry, blockex_1, blockex_2, blockex_3, resadfGeoTransform, blockex_width, blockex_height);//if (GDALRasterIO(merge_band1, GF_Write, 0, 0, blockex_width, blockex_height,blockex_1, blockex_width, blockex_height, GDT_Byte, 0, 0) != CE_None) {//message = "Error : get tmp image band buffer failure, feature id :" + icur;// = false;if (output_dataset)GDALClose(output_dataset);return false;}if (GDALRasterIO(merge_band2, GF_Write, 0, 0, blockex_width, blockex_height,blockex_2, blockex_width, blockex_height, GDT_Byte, 0, 0) != CE_None) {//message = "Error : get tmp image band buffer failure, feature id :" + icur;//ibreak = false;if (output_dataset)GDALClose(output_dataset);return false;}if (GDALRasterIO(merge_band3, GF_Write, 0, 0, blockex_width, blockex_height,blockex_3, blockex_width, blockex_height, GDT_Byte, 0, 0) != CE_None) {//message = "Error : get tmp image band buffer failure, feature id :" + icur;//ibreak = false;if (output_dataset)GDALClose(output_dataset);return false;}CPLFree(blockex_1);CPLFree(blockex_2);CPLFree(blockex_3);//保存tif，删除tifif (output_dataset) {GDALClose(output_dataset);output_dataset = nullptr;}string sss = SavePath + ccc + ".tif";//BOOL bo = DeleteFileA(sss.c_str());QDir dir;dir.remove(QString::fromStdString(sss));int cc = 0;
}void ImageExtract::DrawGeometry(OGRGeometry *geometry, char* blockex_1, char *blockex_2, char *blockex_3, double af[], int width, int height)
{if (geometry == nullptr) return;OGRwkbGeometryType   gt = geometry->getGeometryType();if (gt == wkbPolygon || gt == wkbPolygon25D || gt == wkbPolygonM || gt == wkbPolygonZM){OGRPolygon *polygon = (OGRPolygon*)geometry;OGRLinearRing *ERing = polygon->getExteriorRing();int num_points = ERing->getNumPoints();double* arr_x = new double[num_points], *arr_y = new double[num_points];for (int i = 0; i < ERing->getNumPoints(); i++){OGRPoint *pt = new OGRPoint();ERing->getPoint(i, pt);int r, c;this->LatLon2ColRow(af, pt->getX(), pt->getY(), r, c);arr_x[i] = c; arr_y[i] = r;}std::vector<OGRPoint> block_points;CalGeometry2Pixels(width, height, num_points, arr_x, arr_y, block_points);delete arr_x, arr_x = nullptr;delete arr_y, arr_y = nullptr;DrawGeometryPixels(block_points, blockex_1, blockex_2, blockex_3, width, height);block_points.clear();for (int i = 0; i < polygon->getNumInteriorRings(); i++){OGRLinearRing *InRing = polygon->getInteriorRing(i);num_points = InRing->getNumPoints();arr_x = new double[num_points];arr_y = new double[num_points];for (int j = 0; j < InRing->getNumPoints(); j++){OGRPoint *pt = new OGRPoint();InRing->getPoint(j, pt);int r, c;this->LatLon2ColRow(af, pt->getX(), pt->getY(), r, c);arr_x[j] = c; arr_y[j] = r;}CalGeometry2Pixels(width, width, num_points, arr_x, arr_y, block_points);delete arr_x, arr_x = nullptr;delete arr_y, arr_y = nullptr;DrawGeometryPixels(block_points, blockex_1, blockex_2, blockex_3, width, height);block_points.clear();}}else if (gt == wkbMultiPolygon || gt == wkbMultiPolygon25D || gt == wkbMultiPolygonM || gt == wkbMultiPolygonZM){OGRMultiPolygon *MultPolygon = (OGRMultiPolygon*)geometry;for (int k = 0; k < MultPolygon->getNumGeometries(); k++){OGRGeometry *pGeo = MultPolygon->getGeometryRef(k);if (pGeo->getGeometryType() == wkbPolygon){OGRPolygon *polygon = (OGRPolygon*)pGeo;OGRLinearRing *ERing = polygon->getExteriorRing();int num_points = ERing->getNumPoints();double* arr_x = new double[num_points], *arr_y = new double[num_points];for (int i = 0; i < ERing->getNumPoints(); i++){OGRPoint *pt = new OGRPoint();ERing->getPoint(i, pt);int r, c;this->LatLon2ColRow(af, pt->getX(), pt->getY(), r, c);arr_x[i] = c; arr_y[i] = r;}std::vector<OGRPoint> block_points;CalGeometry2Pixels(width, height, num_points, arr_x, arr_y, block_points);delete arr_x, arr_x = nullptr;delete arr_y, arr_y = nullptr;DrawGeometryPixels(block_points, blockex_1, blockex_2, blockex_3, width, height);block_points.clear();for (int i = 0; i < polygon->getNumInteriorRings(); i++){OGRLinearRing *InRing = polygon->getInteriorRing(i);num_points = InRing->getNumPoints();arr_x = new double[num_points];arr_y = new double[num_points];for (int j = 0; j < InRing->getNumPoints(); j++){OGRPoint *pt = new OGRPoint();InRing->getPoint(j, pt);int r, c;this->LatLon2ColRow(af, pt->getX(), pt->getY(), r, c);arr_x[j] = c; arr_y[j] = r;}CalGeometry2Pixels(width, height, num_points, arr_x, arr_y, block_points);delete arr_x, arr_x = nullptr;delete arr_y, arr_y = nullptr;DrawGeometryPixels(block_points, blockex_1, blockex_2, blockex_3, width, height);block_points.clear();}}}}
}void ImageExtract::CalGeometry2Pixels(int nRasterXSize, int nRasterYSize, int nPartCount, double *padfX, double *padfY, std::vector<OGRPoint>& points)
{if (!nPartCount)return;for (int j = 1; j < nPartCount; j++) {int iX = static_cast<int>(floor(padfX[j - 1]));int iY = static_cast<int>(floor(padfY[j - 1]));const int iX1 = static_cast<int>(floor(padfX[j]));const int iY1 = static_cast<int>(floor(padfY[j]));double dfVariant = 0.0;double dfVariant1 = 0.0;int nDeltaX = std::abs(iX1 - iX);int nDeltaY = std::abs(iY1 - iY);// Step direction depends on line direction.const int nXStep = (iX > iX1) ? -1 : 1;const int nYStep = (iY > iY1) ? -1 : 1;// Determine the line slope.if (nDeltaX >= nDeltaY) {const int nXError = nDeltaY << 1;const int nYError = nXError - (nDeltaX << 1);int nError = nXError - nDeltaX;// == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,// but if it is == 0, dfDeltaVariant is not really used, so any// value is okay.const double dfDeltaVariant =nDeltaX == 0? 0.0: (dfVariant1 - dfVariant) / static_cast<double>(nDeltaX);while (nDeltaX-- >= 0) {if (0 <= iX && iX < nRasterXSize&& 0 <= iY && iY < nRasterYSize)points.push_back(OGRPoint(iX, iY));dfVariant += dfDeltaVariant;iX += nXStep;if (nError > 0) {iY += nYStep;nError += nYError;}else {nError += nXError;}}}else {const int nXError = nDeltaX << 1;const int nYError = nXError - (nDeltaY << 1);int nError = nXError - nDeltaY;// == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,// but if it is == 0, dfDeltaVariant is not really used, so any// value is okay.double dfDeltaVariant =nDeltaY == 0? 0.0: (dfVariant1 - dfVariant) / static_cast<double>(nDeltaY);while (nDeltaY-- >= 0) {if (0 <= iX && iX < nRasterXSize&& 0 <= iY && iY < nRasterYSize)points.push_back(OGRPoint(iX, iY));dfVariant += dfDeltaVariant;iY += nYStep;if (nError > 0) {iX += nXStep;nError += nYError;}else {nError += nXError;}}}}
}

GDAL 根据图版范围对影像进行提取 (C++篇)相关推荐

基于python实现高分二号遥感影像水体提取与水质反演（黑臭水体与水体富营养化）
高分二号遥感影像水体提取与水质反演水体提取函数--NDWI 基于几何约束提取河流生成shp,方便后续裁剪水体水质反演最终结果水体提取函数--NDWI 水体提取函数water.py impor ...
【GlobalMapper精品教程】041：从多波段影像中提取单波段方法
本实验讲解globalmapper中文24软件从多波段遥感影像中提取单波段方法. 参考阅读:[Envi风暴]Envi插件大全:多波段拆分工具的巧妙使用文章目录一.多波段影像介绍二.单波段提取方法 ...
基于多特征的遥感影像阴影提取（python+matlab代码）
来源: 那些年"南方测绘杯"科技论文比赛做的工作先看效果: 实现思路: 由于遮挡物的存在,会使辐射源(太阳)辐射能量不能到达地面某些区域,这些区域就是遥感影像上的阴影区域,它们通 ...
Radiomics：使用先进的特征分析方法从医学影像中提取更多信息
Radiomics:使用先进的特征分析方法从医学影像中提取更多信息关键字图像 Radiomics 肿瘤内部肿瘤异质性摘要实体肿瘤具有空间上和时间上的异质性.这就限制了基于分子检测的侵入性活 ...
python爬虫数据提取_入门Python爬虫——提取数据篇
原标题:入门Python爬虫--提取数据篇作者: 李菲来源:人工智能学习圈前言在提取数据这一环节,爬虫程序会将我们所需要的数据提取出来.在上一篇文章<入门Python爬虫 -- 解析数据 ...
gdal（1）GF-1影像提取NDVI
GF-1号多光谱影像 band3为红波段 band4为近红外波段话不多说,直接上代码 # -*- coding:utf_8 -*- """ 作者:ZHOU 日期:20 ...
数字图像处理：基于NDVI的遥感影像植被提取
文章目录一.实验任务与内容二.实验原理三.编程思路四.实验成果图五.源代码 ImageProcessing.h ImageProcessing.cpp Main.cpp 一.实验任务与内容 ...
提取地图中道路_非机器学习方法·从遥感影像中提取道路
本科低年级曾经做过一个提取道路线的题目. 提供的数据:老师给我们了一幅学校周边地区的影像,包括RGB+Infrared四个波段. 一开始使用的方法:当时尚青涩的我和小伙伴们一起使用KMEANS.阈值分 ...
遥感影像道路提取:Improved Road Connectivity by Joint Learning of Orientation and Segmentation
文章目录一.道路提取难点二.本文贡献三.方法 3.1.方向学习 3.2.连接优化 3.3.堆叠式多分支模块一.道路提取难点基于遥感影像提取道路有以下难点: 云和树的遮挡地形.天气.地理环境 ...

GDAL 根据图版范围对影像进行提取 (C++篇)

GDAL 根据图版范围对影像进行提取 (C++篇)相关推荐

最新文章

热门文章