COPY command psql– loading large LiDAR Point Dataset

Ok, so the INSERT statements were too numerous for inputing the LiDAR point dataset (about a billion points… .)  They kept crashing the postgres daemon.  So, I used copy from a CSV file:

COPY base.cuy_lidar_all FROM 'c:/path/cuy_lidar_ground_veg.insert' WITH CSV

Keep your fingers crossed… .

Limitations to Trigger Based Unique Constraint

A unique constraint implemented as a trigger checking hashed geometry seems like a good idea, that is until I applied it to a multi-10GB dataset.  Not surprisingly, it starts off fast on inserts, and slows down a lot as time goes on.  So, I thought I’d approach duplicates another way, by deleting them once they exist.  So for my table:

base.cuy_contours_2

I have hashed my geometry, and added an index:

UPDATE base.cuy_contours_2 SET hash = MD5(ST_AsBinary(the_geom));
CREATE INDEX cuy_contours_2_hash_idx ON base.cuy_contours_2(hash);

I will follow by deleting the duplicate geometries by checking the hashed geometry:

DELETE
 FROM base.cuy_contours_2
 WHERE gid NOT IN
 (SELECT        MAX(dup.gid)
 FROM        base.cuy_contours_2 as dup
 GROUP BY    dup.hash
 );

But, after a very long time, I get an error:

ERROR:  could not read block 356366 of relation 1663/17185/12118368: Permission
denied

So, it’s a  down week, and most of my maps are cached, so we’ll take the postgres service down, and run it in single user mode:

postgres --single -D c:\postgre_data CM

And rerun the command on one line, per single user mode, with no semicolon:

DELETE FROM base.cuy_contours_2 WHERE gid NOT IN (SELECT MAX(dup.gid) FROM base.cuy_contours_2 as dup GROUP BY dup.hash );

and hope for the best... .
Oh and thanks to :
http://www.postgresonline.com/journal/index.php?/archives/22-Deleting-Duplicate-Records-in-a-Table.html
For help with this one.
				

Rendering a vegetation surface model using PovRay (cont.)

Below is my metacode for iterating through a grid of images to render the vegetation surface model from povray.  The real code will be implemented in BASH and AWK, although this would be a perfect use of Python, if I knew how to use it… .

PostGIS Tables:

cuy_veg_points:        Load veg points into database and make 2D geometry
cuy_terrain_points:        Load terrain points into database and make 2D geometry
cuy_veg_height:        Use nearest neighbor to calculate tree height into veg point database
e.g. http://www.bostongis.com/?content_name=postgis_nearest_neighbor_generic

cuy_tile_centroids:        Generate tile centroids and load into table– contains coordinates and record area name

Process:

Read cuy_tile_centroids
For each centroid record,
     query view extent for cuy_veg_height
     output temporary tree_coords.inc
     write worldfile for image output from povray
     write povray file with camera position at centroid
     render povray file

Rendering a vegetation surface model using PovRay

Early on, I was working on the problem of rendering a virtual forest based on real LiDAR data in order to do things like a detailed, vegetation included, viewshed.  There are other things we can do with this data, including discovering the boundary of the canopy (and thus find canopy gaps) and rendering a vegetation surface model.

Let’s start with a single tree.  If we have a tree, height x, and we have a tree model, height 1, then we can place a tree in our scene height 1x.  If we want to use Povray to render a surface model, we can use a linear pattern applied to our tree to render what parts of the tree are closer and which are further.

#include "transforms.inc"

#declare Camera_Location = <0,80,0>;
#declare Camera_Lookat   = <0,0,0>;
#include "tree1.inc"

camera {orthographic
 location  Camera_Location
 look_at Camera_Lookat
 right 100
 up 100}

// Unioned box is used to normalize the scaling of the pigment
union {
 object { TREE
 scale 80
 }
 box {-1000,1000}
 pigment {
 gradient x color_map {
 [0 color rgb 1]
 [1 color rgb 0]
 }
 scale <vlength(Camera_Location-Camera_Lookat),1,1>
 Reorient_Trans(x, Camera_Lookat-Camera_Location)
 translate Camera_Location
 }
 finish {ambient 1}
}

My tree is the default tree from POV-Tree:  http://web.archive.org/web/20071101052625/propro.ru/go/Wshop/povtree/download.html.

Here’s my ini file:

All_Console=Off
All_File="false"
Antialias=Off
Bits_Per_Color=16
Bounding=On
Bounding_Threshold=10
Continue_Trace=Off
Create_Histogram=Off
Debug_Console=On
Debug_File="false"
Display=On
Draw_Vistas=On
End_Column=600
End_Row=600
Fatal_Console=On
Fatal_File="false"
Height=600
Jitter=Off
Light_Buffer=On
Output_Alpha=Off
Output_File_Type=n
Output_To_File=On
Quality=0
Remove_Bounds=On
Render_Console=On
Render_File="false"
Split_Unions=On
Start_Column=1
Start_Row=1
Statistic_Console=On
Statistic_File="false"
Verbose=On
Vista_Buffer=On
Warning_Console=On
Warning_File="false"
Width=600

You’ll notice I’m outputing as a 16-bit per channel PNG, so I can capture the full dynamic range of values without loosing too much of the precision in the distance values.  Now, just a little algebra, and we have distance from camera calculated, but I’ll leave that to a future post.

Here’s where I’m going– imagine a PostGIS database full of LiDAR vegetation points and ground points.  A nearest neighbor search and a little arithmetic tells us the local height of the vegetation point.  Now we iterate through a grid extracting sections of this dataset, using a simple query and a little AWK scripting to output the include file, the povray file, and render the image of tree height.  The tree height interpolator is then a simple 3D tree shaped stamp that we scale, rotate, and stamp all across our landscape, thus creating a PovRay generated digital surface model of vegetation.

Edit:

Oh, and I must credit the source for my height mapping:

http://news.povray.org/povray.advanced-users/thread/%3Cweb.49c79d64a563150ecb94416e0@news.povray.org%3E/