Alright, so let’s begin again with some LiDAR processing in PostGIS. I really liked this discussion of multipoint conversions of LiDAR and other point datasets, but we won’t address making things more efficient here (yet). Besides, these aren’t just points– we do want to retain their attribution… .
Now, any regular database user will recognize my nubeness reading this post. If you are that person, please comment. I’m looking to learn… .
In the LiDAR dataset I have, there are elevations for all points, ground points, vegetation points, building points etc. I’d like to know about how tall my vegetation is by calculating the difference between the height of the vegetation point and the nearest ground point. So, I need to know what the nearest ground point it, and difference the two heights. Since I am interested in the nearest ground point that is as close to directly beneath my vegetation point as possible, we’ll compute nearest neighbor in two dimensions.
I’m going to, once again, take advantage of the resources at BostonGIS, and their Nearest Neighbor function. Install and use according to their directions… .
So, I’ve downloaded the state LiDAR dataset for an Ohio County, and want to do some processing on it. It does contain some fields in it’s original LAS form we might want to keep, specifically (looking at our statement to create the table:
CREATE TABLE lidar
(
x numeric,
y numeric,
z numeric,
intensity integer,
tnumber integer,
number integer,
class integer,
id integer,
vnum integer,
)
WITH (OIDS=FALSE);
ALTER TABLE lidar OWNER TO postgres;
I convert from las to csv with the las2txt utility (as before):
#!/bin/bash
x=0
for f in $( ls *.las); do
x=`expr $x + 1`
echo $x $f
las2txt --parse xyzinrcpM -sep komma $f $f.csv
done
I want to copy from CSV using local user permissions. It’s a little less efficient to do this than to have the database directly read from a CSV file, but it has the advantage that I don’t have to muck about giving my database super user file level permissions on my OS. With the following steps, I’m doing a sequential test of each step, so while the last step converted all my las files to CSV, I’m still exploring this process with a single file.
So, creating a script file called copy.txt with the following contents
\copy lidar from '/Users/blaa/Documents/GIS/Vinton/OSIP_VIN_LiDAR/s2035470.las.csv' with csv
We load a single block of LiDAR data in:
psql -U postgres -d Vinton -f copy.txt
(note, we could have used a pipe to accomplish this)
Now to add a geometry field, and populate that with the info in our x and y fields (remember, I don’t want these to be 3D points, but 2D points with a z attribute).
SELECT AddGeometryColumn('public','lidar','the_geom',91023236,'POINT', 2);
UPDATE lidar SET the_geom = ST_SetSRID(ST_MakePoint(x,y),91023236);
It’d be good to have a spatial index and primary key… .
CREATE INDEX lidar_the_geom_idx ON lidar USING gist(the_geom);
ALTER TABLE lidar ADD PRIMARY KEY (vnum);
And view in uDig:
OK. Now, I now that class 1 is ground and class 5 is vegetation, so, to simplify future processing, I will split this out into two tables (although subselects or views might work just as well here… ).
CREATE TABLE lidar_veg AS
SELECT * FROM lidar WHERE class = 5;
CREATE TABLE lidar_ground AS
SELECT * FROM lidar WHERE class = 1;
CREATE INDEX lidar_veg_the_geom_idx ON lidar USING gist(the_geom);
CREATE INDEX lidar_ground_the_geom_idx ON lidar USING gist(the_geom);
ALTER TABLE lidar_veg ADD PRIMARY KEY (vnum);
ALTER TABLE lidar_ground ADD PRIMARY KEY (vnum);
The nearest neighbor function is “pgis_fn_nn”, and takes a few options. We’ll try this first with just one point. The vegetation point is in the middle of these 4. There are ground points all around:
SELECT pgis_fn_nn(b.the_geom, 20,1,4,'lidar_ground','true','vnum','the_geom') FROM lidar_veg AS b WHERE b.vnum = 576455;
pgis_fn_nn
------------------
(576457,4.31045)
(1 row)
Now, I need to capture that point, find out it’s z coordinate, and subtract that from the z coordinate of my vegetation point… . Not tonight… .
