问题说明
相关代码下载地址：https://download.csdn.net/download/gainichengyichongfu/18383717?spm=1001.2014.3001.5503
这里贴的是除多边形之外的函数（如点、线等），要丰富一点的可以去上述地址下载~
我要做的
我的目标是：根据由多张图像拼起来的shp文件，生成每张图像参与拼接的二值图像（ 0 表示不参与拼接，255 表示参与拼接）；最后生成一张拼接大图
如（我这里还经过了去云的操作，所以拼接大图里面没有云了，这些细节就不要在意了哈~）
输入数据
我的输入shp文件如下（这是由多张图像拼起来的结果图，每个多边形表示不同影像的参与拼接的部分）：
里面由一个图层、多个多边形组成，其中一个多边形信息如下：
可以看到数据类型是多边形，对应的影像名称 “F:\test_data\CST-testdata\1-GF1-7\DataWithSimuClouds-Nocolor-shape\GF1B_PMS_E112.7_N23.0_20191207_L1A1227736448.tif”
对应的影像图 “F:\test_data\CST-testdata\1-GF1-7\DataWithSimuClouds-Nocolor-shape\GF1B_PMS_E112.7_N23.0_20191207_L1A1227736448.tif” 如下：
我希望可以生成这个影像对应的二值化栅格图，如下：
下面开始讲述如何实现！
方法实现
main主函数
里面使用到了两个函数
（1）shpread(shpName, ValidMaskList, PolygonPoint);
（2）PolygonRasterize(PolygonPoint, ImgList, ValidMaskList, pdriver);
GDALAllRegister();
GDALDriver *pdriver = GetGDALDriverManager()->GetDriverByName("GTiff");//注册驱动
string strpath = "F:\\test_data\\CST-testdata\\1-GF1-7\\DataWithSimuClouds-Nocolor-shape";
string shpName = strpath + "\\result\\project.shp";
vector<string> ValidMaskList;
vector<string> ImgList;
vector<string> PolygonPoint;//云掩膜文件名
shpread(shpName, ValidMaskList, PolygonPoint);
PolygonRasterize(PolygonPoint, ImgList, ValidMaskList, pdriver);
（1）shpread(shpName, ValidMaskList, PolygonPoint);
里面使用到的函数 void shpread(string file_path_name, vector &Masklist, vector &PolygonPoint);--------这个函数的目的是：将多边形的各个拐点读取出来，存在point容器里面；因为后面需要多边形的所有拐点来生成二值的mask栅格文件
shpread函数里面有用到了函数 int Get_Polygon(OGRGeometry *poGeometry, vector& OuterRing, vector<vector>& InteriorRing);，我都放在下面了
声明：由于我确定我的shp文件存储的多个多边形，所以不必要的代码我都注释掉了，不必要的输出也省略了
//读shp文件到对应的point容器
void shpread(string file_path_name, vector<string> &Masklist, vector<string> &PolygonPoint)
{
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); // 支持中文路径
CPLSetConfigOption("SHAPE_ENCODING", ""); //属性表支持中文字段
GDALDataset *poDS = (GDALDataset*)GDALOpenEx(file_path_name.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
if (poDS == NULL)
{
printf("Open failed.\n");
exit(1);
}
//cout << "Open successfully!" << endl;
//获取图层数量
//cout << "<--------获取图层数量-------->" << endl;
//int LayerCount = poDS->GetLayerCount();
//cout << "图层数量: " << LayerCount << endl;
//获取shp图层
//cout << "<--------获取shp图层-------->" << endl;
OGRLayer *poLayer = poDS->GetLayer(0); // 根据序号获取相应shp图层，这里表示第一层
//OGRLayer *poLayer = poDS->GetLayerByName("point"); //根据名称获取相应图层
//获取当前图层的属性表结构
//cout << "<--------获取当前图层的属性表结构-------->" << endl;
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
//重置要素读取顺序
//cout << "<--------重置要素读取顺序-------->" << endl;
poLayer->ResetReading(); // ResetReading() 函数功能为把要素读取顺序重置为从第一个开始
//设置要素指针
//cout << "<--------设置要素指针-------->" << endl;
OGRFeature *poFeature; //用于获取图层上的要素
//创建文件存放数据
string shpname = string(file_path_name).substr(string(file_path_name).find_last_of("/") + 1);
shpname = shpname.substr(0, shpname.find_last_of(".")) + "_";
string attrfile = shpname + "attrfile_shp.txt"; // 属性文件
string datafile = shpname + "datafile_shp.txt"; // 数据文件
//cout << attrfile << ' ' << datafile << endl;
fstream tmpfile(attrfile, ios::out);
fstream tmpdatafile(datafile, ios::out);
int num = 0; //用于标记第几个要素
string finename;
while ((poFeature = poLayer->GetNextFeature()) != NULL)
{
//cout << "图层属性->" << endl;
int aa = poFDefn->GetFieldCount();
for (int iField = 0; iField < poFDefn->GetFieldCount(); iField++)
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
//cout << poFieldDefn->GetNameRef() << ": "; // 属性表字段名
switch (poFieldDefn->GetType()) // 不同数据类型按不同方式，根据字段号获取相应字段的数据
{
case OFTInteger:
//printf("%d,", poFeature->GetFieldAsInteger(iField));
tmpfile << poFeature->GetFieldAsInteger(iField) << "\t"; // 写入文件
break;
case OFTInteger64:
//printf(CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64(iField));
tmpfile << poFeature->GetFieldAsInteger64(iField) << "\t"; // 写入文件
break;
case OFTReal:
//printf("%.3f,", poFeature->GetFieldAsDouble(iField));
tmpfile << poFeature->GetFieldAsDouble(iField) << "\t"; // 写入文件
break;
case OFTString:
//printf("%s,", poFeature->GetFieldAsString(iField));
tmpfile << poFeature->GetFieldAsString(iField) << "\t"; // 写入文件
finename = poFeature->GetFieldAsString(iField);
finename = finename + ".mask";
Masklist.push_back(finename);
break;
default:
//printf("%s,", poFeature->GetFieldAsString(iField));
tmpfile << poFeature->GetFieldAsString(iField) << "\t"; // 写入文件
break;
}
/*printf("\n");*/
}
tmpfile << endl;
string polypoint = "POLYGON ((";
OGRGeometry *poGeometry = poFeature->GetGeometryRef();
//判断当前要素的几何类型，是否为点图层，利用 getGeometryType() 函数获取要素类型
if (poGeometry != NULL)
{
auto GeometryType = wkbFlatten(poGeometry->getGeometryType()); // getGeometryType() 返回的类型可能会有2.5D类型，通过宏 wkbFlatten 转换为2D类型
if (GeometryType == wkbPoint) return;
else if (GeometryType == wkbLineString) return;
else if (GeometryType == wkbMultiPoint) return;
else if (GeometryType == wkbMultiLineString) return;
else if (GeometryType == wkbGeometryCollection) return; //cout << "几何体集合图层要素" << endl;
else if (GeometryType == wkbNone) return; //cout << "该数据只有属性表" << endl;
else if (GeometryType == wkbPolygon) // 3
{
//cout << "多边形图层要素" << endl;
vector<XYZInfo> OuterRing;
vector<vector<XYZInfo>> InteriorRing;
Get_Polygon(poGeometry,OuterRing, InteriorRing);
// 将数据写入文件
tmpdatafile << "<----- 第" + to_string(1 + num++) + "个要素（多边形） ---->" << endl;
tmpdatafile << "外环数据" << endl;
//cout << "数据个数：" << OuterRing.size() << endl;
for (int j = 0; j < OuterRing.size(); j++)
{
//tmpdatafile << OuterRing[j].x << "\t" << OuterRing[j].y << "\t" << OuterRing[j].z << endl; // 写入文件
//tmpdatafile << OuterRing[j].x << "\t" << OuterRing[j].y << "\t" << endl;//不关心Z坐标
polypoint = polypoint + DoubleToString(OuterRing[j].x) + " " + DoubleToString(OuterRing[j].y);
if (j < OuterRing.size() - 1) polypoint = polypoint + ",";
if (j == OuterRing.size() - 1) polypoint = polypoint + "))";
}
tmpdatafile << polypoint.c_str() << endl;//不关心Z坐标
tmpdatafile << "内环数据" << endl;
for (int i = 0; i < InteriorRing.size(); i++)
{
for (int k = 0; k < InteriorRing[i].size(); k++)
tmpdatafile << InteriorRing[i][k].x << "\t" << InteriorRing[i][k].y << "\t" << endl;//不关心Z坐标
//tmpdatafile << InteriorRing[i][k].x << "\t" << InteriorRing[i][k].y << "\t" << InteriorRing[i][k].z << endl; // 写入文件
tmpdatafile << endl;
}
}
//else if (GeometryType == wkbMultiPolygon) // 6
//{
// //cout << "多边形集合图层要素" << endl;
// //获取数据
// vector<vector<XYZInfo>> OuterRingVec;
// vector<vector<vector<XYZInfo>>> InteriorRingVec;
// Get_MultiPolygon(poGeometry, OuterRingVec, InteriorRingVec);
// // 将数据写入文件
// OGRMultiPolygon *MultiPolygon = (OGRMultiPolygon *)poGeometry;
// int NumMPolygon = MultiPolygon->getNumGeometries();
// vector<XYZInfo> xyz;
// tmpdatafile << "<----- 第" + to_string(1 + num++) + "个要素（多边形集合） ---->" << endl;
// for (int k = 0; k < OuterRingVec.size(); k++)
// {
// tmpdatafile << "外环" + to_string(k) << endl;
// for (int j = 0; j < OuterRingVec[k].size(); j++)
// {
// tmpdatafile << OuterRingVec[k][j].x << "\t" << OuterRingVec[k][j].y << "\t" << OuterRingVec[k][j].z << endl; // 写入文件
// }
// for (int i = 0; i < InteriorRingVec[k].size(); i++)
// {
// tmpdatafile << "内环" + to_string(i + 1) << endl;
// for (int j = 0; j < InteriorRingVec[k][i].size(); j++)
// tmpdatafile << InteriorRingVec[k][i][j].x << "\t" << InteriorRingVec[k][i][j].y << "\t" << InteriorRingVec[k][i][j].z << endl; // 写入文件
// }
// tmpdatafile << endl;
// }
//}
}
else
{
printf("no geometry\n");
}
//polypoint = polypoint + "))";
PolygonPoint.push_back(polypoint);
OGRFeature::DestroyFeature(poFeature);
}
tmpfile.close();
GDALClose(poDS);
return;
}
//参数1是外环数据，是一个坐标名到坐标值vec的map，参数2是内环数据，比参数1多了一层编号映射
int Get_Polygon(OGRGeometry *poGeometry,vector<XYZInfo>& OuterRing, vector<vector<XYZInfo>>& InteriorRing)
{
OGRPolygon *poPolygon = (OGRPolygon*)poGeometry;//将几何结构转换成多边形类型
OGRLinearRing *pOGRLinearRing = poPolygon->getExteriorRing();
//获取外环数据
XYZInfo outring;
for (int i = 0; i < pOGRLinearRing->getNumPoints(); i++)
{
outring.x = pOGRLinearRing->getX(i);
outring.y = pOGRLinearRing->getY(i);
outring.z = pOGRLinearRing->getZ(i);
OuterRing.push_back(outring);
}
//获取内环数据（一个外环包含若干个内环）
for (int i = 0; i < poPolygon->getNumInteriorRings(); i++)
{
vector<XYZInfo> Ringvec;
pOGRLinearRing = poPolygon->getInteriorRing(i);
for (int j = 0; j < pOGRLinearRing->getNumPoints(); j++)
{
XYZInfo intrring;
intrring.x = pOGRLinearRing->getX(j);
intrring.y = pOGRLinearRing->getY(j);
intrring.z = pOGRLinearRing->getZ(j);
Ringvec.push_back(intrring);
}
InteriorRing.push_back(Ringvec);
}
return 0;
}
（2）PolygonRasterize(PolygonPoint, ImgList, ValidMaskList, pdriver);
接下来是函数 void PolygonRasterize(vector PolygonPoint, vector ImgList, vector Masklist, GDALDriver *pdriver);//循环处理shp里面的每个多边形文件，主要函数是其同名函数
里面使用到的同名函数 void PolygonRasterize(string Geomrtry, string ImgList, string MaskList, GDALDriver *pdriver);//根据多边形的顶点生成mask
函数我都放在下面了
void PolygonRasterize(vector<string> PolygonPoint, vector<string> ImgList, vector<string> Masklist, GDALDriver *pdriver)//根据多边形的顶点生成mask
{
for (int i = 0; i < ImgList.size(); i++)
{
PolygonRasterize(PolygonPoint[i], ImgList[i], Masklist[i], pdriver);
}
}
void PolygonRasterize(string Geomrtry, string ImgList, string MaskList, GDALDriver *pdriver)//根据多边形的顶点生成mask
{
GDALDataset* pInDataset = (GDALDataset*)GDALOpen(ImgList.c_str(), GA_ReadOnly);
int Width = pInDataset->GetRasterXSize();
int Height = pInDataset->GetRasterYSize();
double dGeoTrans[6] = {};
pInDataset->GetGeoTransform(dGeoTrans);
const char* cProjectionRef = pInDataset->GetProjectionRef();
GDALDataset* pOutDataset = pdriver->Create(MaskList.c_str(), Width, Height, 1, GDT_Byte, NULL);
pOutDataset->SetGeoTransform(dGeoTrans);
pOutDataset->SetProjection(cProjectionRef);
OGRGeometryH hGeometry;
char *pWKT = (char*)Geomrtry.c_str();
OGR_G_CreateFromWkt(& pWKT, NULL, &hGeometry);
int iBand[1] = { 1 };
double dBurnValue[1] = { 255 };
CPLErr err = GDALRasterizeGeometries((GDALDatasetH)pOutDataset, 1, iBand, 1, &hGeometry, NULL, NULL, dBurnValue, NULL, NULL, NULL);
GDALClose(pInDataset);
pInDataset = NULL;
GDALClose(pOutDataset);
pOutDataset = NULL;
pWKT = NULL;
delete[] pWKT;
pWKT = NULL;
}
通过上述可以通过shp的多边形文件，得到对应的二值的栅格文件，如：
若直接裁剪影像也是一样的道理，稍微改改就行了
809
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
