您的位置:首页 > 其它

GDAL矢量栅格化

2010-08-31 22:49 381 查看
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;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签:  null delete float