首页 > 代码库 > las点转为Shapefile文件,获取高程信息

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();//创建ShapefileOGRRegisterAll();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();