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，删除tif
	if (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;
				}
			}
		}
	}
}