Saturday, 11 May 2019

Converting from shapefile or raster formats to LAS?



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?



Answer



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
#include
#include
#include
#include
#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 {
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 \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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...