将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