las点转为Shapefile文件，获取高程信息

将LiDAR点文件转换为Shapefile文件，方便ArcGIS9.3版本操作

const char *pSrcFileName = "D:\\LidarTestData\\1.las";
std::ifstream ifs;ifs.open(pSrcFileName, std::ios::in | std::ios::binary);
if(ifs == NULL)
{
   cout<<"null"<<endl;
}
liblas::ReaderFactory f ;
liblas::Reader reader = f.CreateWithStream(ifs);
liblas::Header const& header = reader.GetHeader();
printf("Points count: %d\n",header.GetPointRecordsCount());
//空间参考
liblas::SpatialReference lasSRef = header.GetSRS();
string sSRS = lasSRef.GetWKT(lasSRef.eCompoundOK,true);
const char *pSRef = sSRS.c_str();
//创建Shapefile
OGRRegisterAll();
const char *shpPath = "D:\\LidarTestData\\test.shp";
OGRDataSource *pODS = NULL;
OGRLayer *pPtLayer = NULL;
OGRSFDriver *pDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");
pODS = pDriver->CreateDataSource(shpPath);
OGRSpatialReference *pSRS = new OGRSpatialReference(pSRef);
pPtLayer = pODS->CreateLayer("Point1",pSRS,wkbPoint);
char * pZFiledName = "ZValue";
OGRFieldDefn pZField(pZFiledName,OFTReal);
pPtLayer->CreateField(&pZField,1);
while(reader.ReadNextPoint())
{
	double x = reader.GetPoint().GetX();
	double y = reader.GetPoint().GetY();
	double z = reader.GetPoint().GetZ();
         OGRPoint pt(x,y);
	OGRFeature *pFeature = OGRFeature::CreateFeature(pPtLayer->GetLayerDefn());
	pFeature->SetGeometry(&pt);
	pFeature->SetField(pZFiledName,z);
	if( pPtLayer->CreateFeature(pFeature) != OGRERR_NONE )
         {
		printf("Failed to create feature in shapefile.\n");
                  exit(1);
         }
	OGRFeature::DestroyFeature(pFeature);
}
pPtLayer->SyncToDisk();

las点转为Shapefile文件，获取高程信息