3

I am looking for a solution to convert either from ESRI shapefile format to LAS or from a common raster format to LAS. I have access to following software:

  • ArcGIS 10.2 for Desktop (arcinfo with spatial analyst and 3d analyst extensions),
  • Global Mapper,
  • LASTools (unlicensed),
  • QGIS 2.8,
  • Whitebox GAT.

I have the functionality in the LASTools shp2las.exe tool, but as I do not have a licence it can only handle 2 million points before it adds XYZ noise. Due to the large volume of data I am working with, it is impractical to cut it into thousands of smaller pieces in order to process with tiles under 2 million points.

There was the potential of a workaround using the text2las tool, but outputs from this tool are also unacceptable due to a stepping/contouring effect. Despite efforts to ensure no rounding or truncation of input data occurs during transformation.

Is there a way through open source or free software to successfully convert between these formats?

Or am I looking at buying a licence for FME or LASTools or some other software?

closed as too broad by PolyGeo Dec 11 '17 at 22:26

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

5

Your features must be points, you can convert raster to point using Raster to Point.

With AddXY tool you can add the coordinates of X, Y (and if present Z) to the points, if not there should be some sort of elevation field already there or you're up the creek...

Start ArcMap, turn the display off (no need to refresh), add the points, format the fields and just keep X,Y,Z, open the attribute table (hint: you can reorder fields by dragging in ArcGis 10+ but they must be in the right order before export). Export the table and use LAStools XYZ2las to convert to a las file.

Apparently Global Mapper will do this with Export Raster and Elevation Data, I do not know if this tool expects Raster or Vector format... there's not much help in the 'help'.

If you are interested here is the (working) C++ code.

#include "stdafx.h"
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <windows.h>
#include <string>
#pragma pack(push,1) // set structure packing to 1

struct LAShdrFirst
{
    char FileSig[4];                 
    unsigned __int16 FileSource;                      
    unsigned __int16 Reserved_Unused;                 
    unsigned __int32 Project_ID_Data1;                 
    unsigned __int16 Project_ID_Data2;                
    unsigned __int16 Project_ID_Data3;               
    char Project_ID_Data4[8];                       
    unsigned char Version_Major;                    
    unsigned char Version_Minor;                    
    char System_ID[32];                             
    char Software[32];                              
    unsigned __int16 FC_Day, FC_Year, Header_Size;    
    unsigned __int32 Offset_to_Data;                    
    unsigned __int32 VarLenRecs;                        
    unsigned char Pt_DataFormat;                   
};
struct LAShdrSecond
{
    unsigned __int16 Pt_DataRecLen;                    
    unsigned __int32 PointCount;                        
    unsigned __int32 Point_by_Return_0;                 
    unsigned __int32 Point_by_Return_1;                 
    unsigned __int32 Point_by_Return_2;                 
    unsigned __int32 Point_by_Return_3;                 
    unsigned __int32 Point_by_Return_4;                 
};
struct LAShdrSpatial
{
    double Xscale;                                  
    double Yscale;                                  
    double Zscale;                                  
    double Xoffset;                                
    double Yoffset;                                
    double Zoffset;                                 
    double MaxX;                                    
    double MinX;                                    
    double MaxY;                                    
    double MinY;                                    
    double MaxZ;                                    
    double MinZ;                                    
};
struct LASrec0
{
    int X; 
    int Y; 
    int Z; 
    unsigned short Intensity; 
    char DataByte; 
    unsigned char Classification; 
    char ScanAngle; 
    char UserData;  
    unsigned short PointSourceID; 
;
#pragma pack (pop) // return packing to whatever it was

const int MaxTxtSize = 256;

void Convert(char* Directory,char* InFileName)
{
    printf("Convert %s\\%s\n",Directory,InFileName);
    char* nName = (char*) malloc(MaxTxtSize);
    for (int i = 0;i<MaxTxtSize;i++)
    {
        if (InFileName[i] == '.')
        {
            nName[i] = '_';
        }
        else if (InFileName[i] == '\0')
        {
            // end of string
            nName[i] = '.';
            nName[i+1] = 'l';
            nName[i+2] = 'a';
            nName[i+3] = 's';
            nName[i+4] = '\0';
            break;
        }
        else
        {
            nName[i] = InFileName[i];
        }
    }
    printf("%s to %s\n",InFileName,nName);
    LAShdrFirst* First = (LAShdrFirst*) calloc(sizeof(LAShdrFirst),1);
    LAShdrSecond* Second = (LAShdrSecond*) calloc(sizeof(LAShdrSecond),1);
    LAShdrSpatial* Spatial = (LAShdrSpatial*) calloc(sizeof(LAShdrSpatial),1);
    LASrec0* Record = (LASrec0*) calloc(sizeof(LASrec0),1);
    char* Empty4 = (char*) calloc(4,1);

    First->FileSig[0]='L';
    First->FileSig[1]='A';
    First->FileSig[2]='S';
    First->FileSig[3]='F';

    First->Offset_to_Data=229;
    First->Pt_DataFormat = 0;
    First->Version_Major = 1;
    First->Version_Minor = 2;
    Second->Pt_DataRecLen=20;
    Spatial->Xscale = 0.01;
    Spatial->Yscale = 0.01;
    Spatial->Zscale = 0.01;

    char* InFilePath = (char*) malloc(MaxTxtSize);
    char* OutFilePath = (char*) malloc(MaxTxtSize);
    std::string FileLine;
    char* FileData;
    sprintf(InFilePath,"%s\\%s",Directory,InFileName);
    sprintf(OutFilePath,"%s\\%s",Directory,nName);
    double pX,pY,pZ;
    unsigned __int32 PointCount = 0;
    unsigned __int16 Intensity  = 0;
    char* TextValues = (char*) malloc(32);
    bool ExtentStarted = false;
    for (int i = 0;i < 32;i++)
    {
        TextValues[i] = '\0';
    }
    int vPnt = 0;
    int MaxLen;
    int pCnt,TVp;

    std::ifstream InFile (InFilePath,std::ios::in);
    std::ofstream OutLasFile(OutFilePath,std::ios::out | std::ios::binary);
    OutLasFile.write((char*)First,sizeof(LAShdrFirst));
    OutLasFile.seekp(First->Offset_to_Data,std::ios::beg); // seek to first record

    do
    {
        getline(InFile,FileLine);
        MaxLen = FileLine.size();
        if (MaxLen > 0)
        {
            FileData = new char[MaxLen+1];
            std::copy(FileLine.begin(),FileLine.end(),FileData);
            FileData[MaxLen] = '\0'; // null terminate VERY IMPORTANT

            vPnt = 0;
            TVp = 0;
            for (pCnt = 0;pCnt <= MaxLen;pCnt++)
            {
                if (FileData[pCnt] == ' ' || FileData[pCnt] == ',' || FileData[pCnt] == '\t' || FileData[pCnt] == '\0') // comma or end string.
                {
                    if (vPnt == 0)
                        pX = atof(TextValues);
                    else if (vPnt == 1)
                        pY = atof(TextValues);
                    else if (vPnt == 2)
                        pZ = atof(TextValues);
                    else if (vPnt == 3)
                        Intensity = atoi(TextValues);
                    TVp = 0;
                    vPnt++;
                    for (int i = 0;i < 32;i++)
                    {
                        TextValues[i] = '\0';
                    }
                }
                else
                {
                    TextValues[TVp] = FileData[pCnt];
                    TVp++;
                }
            }

            if (ExtentStarted)
            {
                Spatial->MinX =(pX < Spatial->MinX) ? pX : Spatial->MinX;
                Spatial->MaxX =(pX > Spatial->MaxX) ? pX : Spatial->MaxX;
                Spatial->MinY =(pY < Spatial->MinY) ? pY : Spatial->MinY;
                Spatial->MaxY =(pY > Spatial->MaxY) ? pY : Spatial->MaxY;
                Spatial->MinZ =(pZ < Spatial->MinZ) ? pZ : Spatial->MinZ;
                Spatial->MaxZ =(pZ > Spatial->MaxZ) ? pZ : Spatial->MaxZ;
            }
            else
            {
                Spatial->MinX = pX;
                Spatial->MaxX = pX;
                Spatial->MinY = pY;
                Spatial->MaxY = pY;
                Spatial->MinZ = pZ;
                Spatial->MaxZ = pZ;
                ExtentStarted = true;
            }
            Record->X = pX * 100;
            Record->Y = pY * 100;
            Record->Z = pZ * 100;
            Record->Intensity = Intensity;
            OutLasFile.write((char*)&Record->X,4);
            OutLasFile.write((char*)&Record->Y,4);
            OutLasFile.write((char*)&Record->Z,4);
            OutLasFile.write((char*)&Record->Intensity,2);
            OutLasFile.write((char*)&Record->DataByte,1);
            OutLasFile.write((char*)&Record->Classification,1);
            OutLasFile.write(Empty4,4);
            PointCount++;
            delete[] FileData;
        }
    } while (!InFile.eof());
    OutLasFile.flush();

    Second->Point_by_Return_0 = PointCount;
    Second->PointCount = PointCount;

    OutLasFile.seekp(105,std::ios::beg);
    OutLasFile.write((char*)&Second->Pt_DataRecLen,2);
    OutLasFile.write((char*)&Second->PointCount,4);
    OutLasFile.write((char*)&Second->Point_by_Return_0,4);
    OutLasFile.write(Empty4,4);// not writing points-by-return
    OutLasFile.write(Empty4,4);
    OutLasFile.write(Empty4,4);
    OutLasFile.write(Empty4,4);

    OutLasFile.write((char*)Spatial,sizeof(LAShdrSpatial));

    InFile.close();
    OutLasFile.close();
    free(TextValues);
    free(nName);
    free(InFilePath);
    free(OutFilePath);
    free(First);
    free(Second);
    free(Spatial);
    free(Record);
    free(Empty4);
}

int _tmain(int argc, _TCHAR* argv[])
{
    if (argc == 1)
    {
        printf("XYZI comma or space delimeted format to LAS converter version 1.1.\n");
        printf("\n");
        printf("Tool written in C++ by Michael Miles-Stimson of RPS Spatial\n");
        printf("\n");
        printf("Converts X,Y,Z(,I) formatted file into LAS version 1.2 format.\n");
        printf("Please ensure there is no header line in the XYZI file\n");
        printf("Before running this tool.\n");
        printf("\n");
        printf("This tool converts every file in the input folder that has the\n");
        printf("extension .XYZI or .XYZ, naming the file with the full name appending\n");
        printf(".las, i.e. abcd.xyzi becomes abcd_xyzi.las.\n");
        printf("\nAll points will have classification of 0\n");
        printf("\n");
        printf("Usage: XYZI_To_LAS <Input Folder>\n");
        return 0;
    }
    char* fName = (char*) malloc(MaxTxtSize);
    char* dName = (char*) malloc(MaxTxtSize);
    wcstombs(dName,argv[1],MaxTxtSize);

    wchar_t* wtf = new wchar_t[MaxTxtSize];
    swprintf(wtf,MaxTxtSize,L"%s\\*.XYZI",argv[1]);
    WIN32_FIND_DATA fData;
    HANDLE hFind;
    hFind = FindFirstFile(wtf,&fData);
    if (hFind != INVALID_HANDLE_VALUE)
    {
        wcstombs(fName,fData.cFileName,MaxTxtSize);
        //printf("First File %s\n",fName);
        Convert(dName,fName);

        while (FindNextFile(hFind,&fData) != 0)
        {
            wcstombs(fName,fData.cFileName,MaxTxtSize);
            Convert(dName,fName);
            //printf("> %s\n",fName);
        }
    }

    swprintf(wtf,MaxTxtSize,L"%s\\*.XYZ",argv[1]);
    hFind = FindFirstFile(wtf,&fData);
    if (hFind != INVALID_HANDLE_VALUE)
    {
        wcstombs(fName,fData.cFileName,MaxTxtSize);
        Convert(dName,fName);
        while (FindNextFile(hFind,&fData) != 0)
        {
            wcstombs(fName,fData.cFileName,MaxTxtSize);
            Convert(dName,fName);
            //printf("> %s\n",fName);
        }
    }
    delete[] wtf;
    free(fName);
    free(dName);
    return 0;
}

This converts a folder of XYZI files into ASPRS LAS 1.2 with point record type 0. The extension must be XYZI but can be comma, space or tab delimited. There's not much commenting and it's quite sloppy but it was only my second C++ program. Copy and paste the code into a new Visual C++ project (console application), build and then run on the command line.

Note: I have physically embedded only a part of my LASheader.h (2266 lines in this header) into this code so it's possible there's a struct missing, I've given it a quick check and can't see anything obvious.

  • Thank you for the suggestion. I have tried using the XYZtolas / text2las tool before as a work-around, but I have found that it leaves an undesirable "countour-like effect" or "steps" on the output elevation data. – Jordan Walker Mar 19 '15 at 22:45
  • The steps are likely to come from your input data, that or your elevations are integer/rounded. Would you be interested in a python solution to export your XYZ values to perhaps avoid rounding issues? – Michael Stimson Mar 19 '15 at 22:47
  • I would have thought that would be the case too, but there is no noise in the input data, and I used a number of different routes to obtain the XYZ text file, none of them rounded decimal places or performed an int function... but the steps are caused as a function of the text2las tool as far as I can tell. I am open to any solution, really. – Jordan Walker Mar 19 '15 at 22:49
  • I wrote my own code in C++ to do this; although we have downloaded LAStools the license that was purchased didn't include me so I can only 'evaluate', hence I wrote my own. According to resources.arcgis.com/en/help/main/10.1/index.html#//… the XYZ is a valid LAS format according to ArcGis. – Michael Stimson Mar 19 '15 at 23:17
  • Is the C++ code something you would be willing to share? If so, how do you usually implement it, and can I batch process with it? If not, I appreciate your help all the same. – Jordan Walker Mar 20 '15 at 0:50
1

The PDAL GDAL reader will read any GDAL-readable raster source and turn it into a point cloud. You can then use other PDAL facilities as necessary filter and process the data.

Not the answer you're looking for? Browse other questions tagged or ask your own question.