C# GDAL基础使用之矢量栅格化处理6

最新推荐文章于 2023-03-15 11:02:43 发布

Prince_Treasure

最新推荐文章于 2023-03-15 11:02:43 发布

阅读量855

点赞数 1

分类专栏： GDAL

本文链接：https://blog.csdn.net/prince999999/article/details/118668726

版权

GDAL 专栏收录该内容

31 篇文章 9 订阅

订阅专栏

该代码段展示了使用GDAL和OGR库将矢量数据转换为栅格的过程。它首先注册GDAL和OGR，然后打开输入的矢量文件，获取其地理信息和范围。接着，创建一个指定尺寸的GeoTIFF栅格文件，并设置地理变换。如果矢量数据无地理信息，代码会自定义地理变换。最后，使用GDALRasterizeLayers进行矢量数据的栅格化操作，并保存结果。

摘要由CSDN通过智能技术生成


CString COGRVToR::Run()
{
	GDALAllRegister();
	OGRRegisterAll();
	CPLSetErrorHandler(CPLLoggingHandler);//在刚开始的GDAL调用时
	CString strErr;
	OGRDataSource  *poDS;		
	poDS = OGRSFDriverRegistrar::Open(m_strInputV, FALSE );
	if( poDS == NULL )
	{
		strErr.Format("Open %s failed./n",m_strInputV);
		return strErr;
		//	exit( 1 );
	}		
	OGRLayer* poLayer;	
	poLayer = poDS->GetLayer(0);	
	const char * pszOutFileName =m_strOutputR;
	char * pszProjection;
	OGRSpatialReference * poSRS = poLayer->GetSpatialRef();
	OGREnvelope  psExent;
	poLayer->GetExtent(&psExent);
	int m_nImageWidth =1024;// podataset->GetRasterXSize();
	int m_nImageHeight =1024;// podataset->GetRasterYSize();	
	if ( !poSRS )
	{
		strErr.Format("无法打开%s地理信息",m_strInputV);
		m_nImageWidth = psExent.MaxX;
		m_nImageHeight = psExent.MaxY;
		CString str;
		str.Format("%d ,%d" , m_nImageWidth,m_nImageHeight);
		str = str+strErr;
		AfxMessageBox(str);
		//	AfxMessageBox(strErr);
		//	exit(1);
		//	return strErr;
	}
	else
	{
		poSRS->exportToWkt( &pszProjection );

	}

	
	const char *pszFormat = "GTiff";
	GDALDriver *poDriver;
	poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);		
//	if( poDriver == NULL )	
//		exit( 1 );	
	GDALDataset *poNewDS = poDriver->Create( pszOutFileName, m_nImageWidth,m_nImageHeight,1, GDT_Float32 , NULL );	
	if ( !poSRS )
	{
	//	strErr.Format("因为无地理信息,%s无法栅格化",m_strInputV);
		
	//	AfxMessageBox(strErr);
	//	return strErr;
		//fot jy
		double adfGeoTransform[6];
		adfGeoTransform[0] = psExent.MinX;
		adfGeoTransform[1] = (psExent.MaxX - psExent.MinX)/m_nImageWidth;
		adfGeoTransform[2] = 0;
		adfGeoTransform[3] = psExent.MaxY;
 		adfGeoTransform[4] = 0;
 		adfGeoTransform[5] = (psExent.MinY - psExent.MaxY)/m_nImageHeight;				
		GDALSetGeoTransform( poNewDS, adfGeoTransform );
// 		OGRSpatialReference * poOGRSR;
// 		poOGRSR->SetUTM(50);
// 		poSRS->exportToWkt(&pszProjection);
//  		poNewDS->SetProjection(pszProjection);	
		//for jy

	}
	else
	{
		double adfGeoTransform[6];
		adfGeoTransform[0] = psExent.MinX;
		adfGeoTransform[1] = (psExent.MaxX - psExent.MinX)/m_nImageWidth;
		adfGeoTransform[2] = 0;
		adfGeoTransform[3] = psExent.MaxY;
		adfGeoTransform[4] = 0;
		adfGeoTransform[5] = (psExent.MinY - psExent.MaxY)/m_nImageHeight;				
		GDALSetGeoTransform( poNewDS, adfGeoTransform );
		poNewDS->SetProjection(pszProjection);	
		
	} 	
//	GDALDatasetH poDstDs;
//	poDstDs = (GDALDatasetH)poNewDS;//GDALOpen(pszOutFileName,GA_Update);	
	int * pnbandlist;
	pnbandlist = new int [1];
	pnbandlist[0]=1;
	double *dburnValues;
	dburnValues = new double[1];
	dburnValues[0]=255;
	OGRLayerH * player;
	player = new OGRLayerH[1];
	player[0] = (OGRLayerH)poLayer;
	//	poLayer->GetInfo()
	char **papszOptions = NULL;
	//	papszOptions = CSLSetNameValue( papszOptions, "ALL_TOUCHED", "TRUE" );
	//	papszOptions = CSLSetNameValue( papszOptions, "ATTRIBUTE", "Band1Mean" );	
	papszOptions = CSLSetNameValue( papszOptions, "CHUNKSIZE", "1" );
	papszOptions = CSLSetNameValue( papszOptions, "ATTRIBUTE", m_strPara );	
				
//	AfxMessageBox("start");
//	CPLErr err= GDALRasterizeLayers(poDstDs,1,pnbandlist,1,player,NULL,NULL,dburnValues,papszOptions,NULL,NULL);
 //  void * pTransformArg;
 //  pTransformArg = GDALCreateApproxTransformer();
	void * pTransformArg=NULL;
	void * m_hGenTransformArg=NULL;
	 m_hGenTransformArg = GDALCreateGenImgProjTransformer( NULL,
                                         pszProjection,
                                         (GDALDatasetH)poNewDS,
                                         poNewDS->GetProjectionRef(),
                                         FALSE, 1000.0, 3 );
	pTransformArg = GDALCreateApproxTransformer( GDALGenImgProjTransform,m_hGenTransformArg, 0.125 );

	CPLErr err= GDALRasterizeLayers((GDALDatasetH)poNewDS,1,pnbandlist,1,player,GDALGenImgProjTransform,m_hGenTransformArg,NULL,papszOptions,GDALTermProgresscjs,"矢量栅格化");
//	AfxMessageBox("end");
				//	CPLErr err= GDALRasterizeLayers(po
//	DstDs,1,pnbandlist,1,player,NULL,NULL,NULL,papszOptions,NULL,NULL);
	 GDALDestroyGenImgProjTransformer(m_hGenTransformArg);
	 GDALDestroyApproxTransformer(pTransformArg);
	 GDALClose(poNewDS);   
	OGR_DS_Destroy( poDS );
//	GDALClose( poDstDs );
	delete[]player;
	delete[]pnbandlist;
	delete[]dburnValues;
	strErr = "栅格化成功";
	return strErr;
}