GDAL矢量转栅格

gdal版本1.8.0.

首先尝试使用gdal_rasterize小工具，但是不支持新生成输出栅格文件。官方说1.8版本后已经支持。但仍然无效。

栅格的元数据赋值和仿射变换以后需要进一步研究。

代码实现：

#include "stdafx.h"
#include <iostream>
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include "gdal_alg.h"

using namespace std;

bool run_rasterize(const char* shp,const char* rasterFile);
int _tmain(int argc, _TCHAR* argv[])
{
	if (run_rasterize("..\\gis_data\\bountl.SHP","..\\gis_data\\bountl3.tif"))
	{
		cout<<"ok...\n";
	}
	else
	{
		cout<<"fail...\n";
	}
	
	cin.get();
	return 0;
}

bool run_rasterize(const char* shp,const char* rasterFile)
{
	GDALAllRegister();  
	OGRRegisterAll(); 

	//打开矢量图层
	OGRDataSource *pOgrSrc = NULL;
	pOgrSrc = OGRSFDriverRegistrar::Open( shp, false);

	if (pOgrSrc == NULL)
	{
		return false;
	}
	OGRLayer *pOgrLyr;
	pOgrLyr = pOgrSrc->GetLayer(0);
	
	//
	OGREnvelope env;
	pOgrLyr->GetExtent(&env);

	int m_nImageWidth =1024; 
	int m_nImageHeight =1024;

	OGRSpatialReference *pOgrSRS = NULL;
	pOgrSRS = pOgrLyr->GetSpatialRef();

	char *pPrj = NULL;
	if (pOgrSRS == NULL)
	{
		cout<<"无投影信息...\n";
		m_nImageHeight = (int)env.MaxX;
		m_nImageWidth = (int)env.MaxY;
	}else
	{
		pOgrSRS->exportToWkt(&pPrj);
	}

	const char *pszFormat = "GTiff"; 
	GDALDriver *poDriver = NULL;  
	poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); 

	GDALDataset *poNewDS = poDriver->Create(rasterFile, 
			                                 m_nImageWidth,
											 m_nImageHeight,
											 1, 
											 GDT_Float32 , 
											 NULL ); 
	double adfGeoTransform[6];  
	adfGeoTransform[0] = env.MinX;  
	adfGeoTransform[1] = (env.MaxX - env.MinX)/m_nImageWidth;  
	adfGeoTransform[2] = 0;  
	adfGeoTransform[3] = env.MaxY;  
	adfGeoTransform[4] = 0;  
	adfGeoTransform[5] = (env.MinY - env.MaxY)/m_nImageHeight;                
	GDALSetGeoTransform( poNewDS, adfGeoTransform );

	if (pOgrSRS != NULL)
	{
		poNewDS->SetProjection(pPrj);
	} 
 
	int * pnbandlist = NULL;  
	pnbandlist = new int [1];  
	pnbandlist[0]=1;  
	double *dburnValues = NULL;  
	dburnValues = new double[3];  
	dburnValues[0]=255;  
	dburnValues[1]=111;  
	dburnValues[2]=34;  

	OGRLayerH * player;  
	player = new OGRLayerH[1];  
	player[0] = (OGRLayerH)pOgrLyr;  
	 
	char **papszOptions = NULL;  
 	papszOptions = CSLSetNameValue( papszOptions, "CHUNKSIZE", "1" );  
	papszOptions = CSLSetNameValue( papszOptions, "ATTRIBUTE", "CODE" );   
	
	void * pTransformArg=NULL;  
	void * m_hGenTransformArg=NULL;  
	m_hGenTransformArg = GDALCreateGenImgProjTransformer( NULL,  
		pPrj,  
		(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,
		dburnValues,
		papszOptions,
		GDALTermProgress,
		"vector2raster");
	
	GDALDestroyGenImgProjTransformer(m_hGenTransformArg);  
	GDALDestroyApproxTransformer(pTransformArg);  
	GDALClose(poNewDS);     
	 
	OGRDataSource::DestroyDataSource(pOgrSrc);
	delete[]player;  
	delete[]pnbandlist;  
	delete[]dburnValues; 

	return true;
}

效果图：

图1. 线要素矢量图

图2. 栅格化后图