将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文件,获取高程信息,布布扣,bubuko.com
原文:http://blog.csdn.net/rybgis/article/details/26017889