For a variety of reasons, I want to load my whole LiDAR vegetation dataset into PostGIS.  I have it in raw, unclassified LAS files, and broken into classified, space delimited text files as well.  The unclassified LAS is nice, if I ever want to hone the data, but for now, we’ll assume the classified data are correctly classified by someone smarter than I am (a likely scenario).

First I wanted to concatenate the text files (here in Windows Command):

copy *.xyz cuy_lidar_veg.txt

Now, I really wanted them in CSV form to make my insert statements easy, so I cheated and used lastools to do the job, assuming it would be more efficient than a stream based script, like awk or sed:

txt2las -i cuy_lidar_veg.txt -o cuy_lidar_veg.las -parse xyz
las2txt -i cuy_lidar_veg.las -o cuy_lidar_veg.csv -sep comma-parse xyz

converting from this:

2114483.860 612559.340 794.720
2114470.510 612558.420 794.420
2114449.970 612565.010 793.900
2114451.540 612569.640 793.800
2114481.010 612571.740 794.390
2114473.230 612572.950 794.130 ...

to this:

2114483.86,612559.34,794.72
2114470.51,612558.42,794.42
2114449.97,612565.01,793.9
2114451.54,612569.64,793.8
2114481.01,612571.74,794.39
2114473.23,612572.95,794.13 ...

I need a table to accept the data so in psql:

CREATE TABLE base.cuy_veg_lidar(
 x numeric,
 y numeric,
 z numeric
);

Now to use awk to format a series of insert statements and pass them through psql right into the database:

less cuy_lidar_veg.csv | awk '{ print "INSERT INTO base.cuy_veg_lidar (x,y,z) VALUES (" $0 ");" };' | psql -U postgres -d CM

This creates a series of command which look like this:

INSERT INTO base.cuy_veg_lidar (x,y,z) VALUES (2114120.67,612570.62,791.11);

Oh, and this last step I ran in CYGWIN command prompt so I could easily access my favorite unix utilities… .  I’m still a nube in the Windows world.  I need my crutches… .

Why so much work when I could just use Postgres’ COPY command?  Permissions, permissions, permissions… .  I like my database to have very little file access ability… .

Next we add a geometry column for 3 dimensional points:

SELECT AddGeometryColumn('base','cuy_veg_lidar','the_geom',9102722,'POINT', 3);

And populate that with our XYZ information:

update base.cuy_veg_lidar set the_geom = ST_SetSRID(ST_MakePoint(x,y,z),102722);