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;
}
}
}
}
}
扫一扫
评论
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
查看更多评论
添加红包
