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.

Loading LiDAR data in PostGIS

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);

Unique constraint on geometry continued

An actual example of using the hashed geometry to uniquely constrain geometry (nothing really new here, just an example of using what we learned in the last post with a new problem):

I’ve got a contours dataset from the county where I work that is in 20 different shapefiles ranging in size from 500MB to 2GB.  I want to put them into a PostGIS database in a single table, so that I can use PostGIS to slice and dice them into reasonable chunks for use by our engineers in AutoDesk products, by ArcGIS users in several divisions, and to improve spatial indexing when dump them into one table in our spatial database as well (ala Regina Obe/Boston GIS’s Map Dicing).  This is a really detailed LiDAR-based contour dataset which used breaklines to make a really detailed and accurate product– which sadly is unusable in its current form because it takes too long to load any given section.

But, today, no slicing and dicing just yet.  I loaded all the data and discovered duplicate geometries.  It seems that some of these blocks that the data come in have overlaps.  I don’t have a screen shot of it, but when I viewed the data in uDig, you can see the duplicate geometries because the transparency of the contour lines in the duplicated areas is less transparent, i.e. the lines are darker against a white background.  I could search for and remove duplicates, but that seems heavy handed.  I’m going to recreate the table with a unique constraint applied to the hashed geometry (and hope for no hash collisions).

First we create the table:

CREATE TABLE base.cuy_contours_2
(
 gid serial NOT NULL,
 elevation smallint,
 shape_leng real,
 the_geom geometry,
 hash character(32),
 CONSTRAINT cuy_contours_2_pkey PRIMARY KEY (gid),
 CONSTRAINT geom_uniq UNIQUE (hash),
 CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 4),
 CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'LINESTRING'::text OR the_geom IS NULL),
 CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 9102722)
)
WITH (OIDS=FALSE);
ALTER TABLE base.cuy_contours_2 OWNER TO postgres;

CREATE INDEX cuy_contours_2_the_geom_idx
 ON base.cuy_contours_2
 USING gist
 (the_geom);

Now we create our trigger that creates a hashed version of the geometry, “the_geom”, in a column called “hash”.

DROP TRIGGER hasher on base.cuy_contours_2;
CREATE OR REPLACE FUNCTION hash_geom() RETURNS TRIGGER AS $hasher$
BEGIN
 IF(TG_OP='INSERT') THEN

 UPDATE base.cuy_contours_2 SET hash = MD5(ST_AsBinary(the_geom));

 END IF;
 RETURN NEW;
END;
$hasher$ LANGUAGE plpgsql;

CREATE TRIGGER hasher AFTER INSERT ON base.cuy_contours_2 FOR EACH STATEMENT EXECUTE PROCEDURE hash_geom();

And finally we add a unique constraint:

ALTER TABLE base.cuy_contours_2 ADD CONSTRAINT geom_uniq UNIQUE (hash);

Now, sadly, I have to use windows where I work, so the following is Windows Command language to load all my shapefiles into the PostGIS database.  That said, I think the command, in this rare case, is more elegant than similar code in, say, BASH scripting:

for %f in (*.shp) do shp2pgsql -s 102722 -a -S -N skip %f base.cuy_contours_2 | psql -h bobs_server -d CM -U postgres

using the shp2pgsql -a flag for append, -S to enforce simple geometries, -N skip to ensure that we keep all records with non-null geometries to compensate for potential errors in the input dataset).

PostGIS Gripe—Limits to Postgre’s B-tree indexing—Followup

In an earlier post, I griped about the limits to Postgres B-tree indexing, in that for large PostGIS geometries, I couldn’t create a unique constraint for PostGIS geometry.  Abe Gillespie (blame me not him for any mistakes  in what follows), suggested that I hash the geometry and create an index on that.  In a follow-up e-mail, he suggested a series of tests in a trigger statement, first test on the bbox, then on the hashed geometry, then finally on the geometry itself upon insert, to determine if the geometry is unique or not.

Well, I’ve completed the first suggestion, but not yet the second.  Here’s my trigger statement:

DROP TRIGGER dissolver on public.reservation_boundaries_public_private_cm;
CREATE OR REPLACE FUNCTION dissolve_bounds() RETURNS TRIGGER AS $dissolver$
BEGIN
   IF(TG_OP='INSERT') THEN

     UPDATE reservation_boundaries_public_private_cm SET hash = MD5(ST_AsBinary(the_geom));

   END IF;
 RETURN NEW;
END;
$dissolver$ LANGUAGE plpgsql;

CREATE TRIGGER dissolver BEFORE AFTER INSERT ON public.reservation_boundaries_public_private_cm FOR EACH STATEMENT EXECUTE PROCEDURE dissolve_bounds();

So, each time I insert a new geometry, the field “hash” is updated with an MD5 hashed version of the geometry.  Now we want to test against this so we add a constraint:

ALTER TABLE public.reservation_boundaries_public_private_cm ADD CONSTRAINT geom_uniq UNIQUE (hash);

And it works a charm. The only problem is that if I were to get a hash collision, I would reject a record and not mean to. I’m not sure the statistical chances of this, but implementing the above flow chart would prevent such a possibility.  If my reading is right, the chances of a hash collision is nearly 1 in 2^128, so maybe I won’t worry this year… .

PostGIS Triggers

Yup, I’m going for hits with the above title.  So, I’ve been playing with triggers in PostGIS (PostgreSQL) trying to streamline a process that I go through quarterly.  We have a geographic boundary layer for our public entity:

The green area is the area with public access.  The purple without.  I want to maintain a version of the boundary with this very simple distinction, and another version of it in which the boundary is dissolved.  Triggers seem ideal here.  When I update the first, I get a version of the second (see below)

The command to do this is as follows:

CREATE TABLE public.reservation_boundaries_public_private_cm_dissolved
 AS
 SELECT res, ST_Union(the_geom) AS the_geom
  FROM public.reservation_boundaries_public_private_cm
   GROUP BY res;

So far, so good (sorry about the awful table names).  Now it needs to run each time there is an insert on the original table.  Here is what I came up with:

DROP TRIGGER dissolver ON public.reservation_boundaries_public_private_cm;

CREATE OR REPLACE FUNCTION dissolve_bounds() RETURNS TRIGGER AS $dissolver$

BEGIN
 IF(TG_OP='INSERT') THEN

 DROP TABLE IF EXISTS public.reservation_boundaries_public_private_cm_dissolved;

 CREATE TABLE public.reservation_boundaries_public_private_cm_dissolved
 AS
 SELECT res, ST_Union(the_geom)
   AS the_geom
   FROM public.reservation_boundaries_public_private_cm
   GROUP BY res;

 END IF;
 RETURN NEW;
END;

$dissolver$ LANGUAGE plpgsql;

CREATE TRIGGER dissolver
  AFTER INSERT ON public.reservation_boundaries_public_private_cm
  FOR EACH STATEMENT EXECUTE PROCEDURE dissolve_bounds();

My only complaint here is speed– it takes about 45 seconds (instead of close to 0.3 seconds) now to insert all my records.  I hoped that I could speed it up with the addition of

FOR EACH STATEMENT EXECUTE PROCEDURE

instead of

FOR EACH RECORD EXECUTE PROCEDURE

but the difference in speed is undetectable.  Any suggestions on how to speed this puppy up?

No Image– a very simple solution for calculating viewshed

Not much to say here. I had a bit of inspiration, remembering a Christoph Hormann tutorial http://www.imagico.de/pov/icons.html on making transparent icons in povray and realized that what he was describing was the final piece in the povray for viewsheds milleau: no_image option for an object forces it to interact with the scene in every way, but to be seen by the camera. I’ll post the revised code coming up, but it’s pretty simple. Add the following option inside the code to the tree object:
no_image .

http://www.povray.org/documentation/view/3.6.1/328/

Another Blog

Just as I’ve started getting good traffic on this blog, a collegue and I have started another one called Anantasesha, largely on the GIS stack we’re putting together for internal use, but also may include some interesting other bits of code and musings. Check it out:
http://anantasesha.wordpress.com

POV-ray for GIS Visualization

This is a “wouldn’t it be nice” post– i.e. a I-don’t-have-time-to-write-it-and-everything-else-on-my-list post. Wouldn’t it be nice if in addition to returning GeoRSS, WFS, OGC, KML, etc. GeoServer returned a POV-Ray file? It could be returned as part of a simple GET post like this:

http://localhost:8080/geoserver/wms?bbox=2100506,575982,2290532,739026&styles=&Format=application/povray&request=GetMap&version=1.1.1&layers=base:contours&width=800&height=644&srs=EPSG:102722&camera=2195519,657504,1000&look_at=2195029,657004,500

Vector symbolization could be a POV-Ray approximation of the WFS render styles using coordinates and POV primitives. Also, a PHP script could consume this return, and run a command line version of POV to render the file in a webserver implementation. This second part I’ll probably actually write with just a DEM and image overlay, but it would be cool to have the vector symbolization as well as part of another GeoServer protocol.  Of course, I recognize this is not a OGC protocol… .  Alternatively, a PHP script to consume WKT and convert it to a POV-Ray file might be easier than incorporating it in GeoServer.

PostGIS Gripe—Limits to Postgre’s B-tree indexing

Is this a gripe about PostGIS or PostgreSQL?  Hard to say.  At my office, we are migrating to a PostGIS—GeoServer—Client-of-your-choice stack.  As usual, I have to try the hard stuff first.  So, rather than keeping the 2ft contour lines for our 7 county area in shapefiles and doing some indexing on them there, I’m attempting to load them into PostGIS all in a single table.  Make any sense?  Maybe not.  Honestly, I really don’t know.  In doing so though, PostGIS has done pretty well.  I loaded all 9.6 GB of 2ft contours from one county into the database, served it up directly in uDig, and then through GeoServer to OpenLayers and Google Earth.  For viewing the whole county, it was slow, but I’d never let a user view the 2ft contour dataset all in one go.  Zoomed in, however, it did pretty well, in fact, I’d say considering there was no caching turned on for OpenLayers, it did splendidly.

First, how I loaded it:

I used shp2pgsql.  My contours were in a bunch of little chunks derived from a state DEM dataset, so I iterated through:

Create the table:

shp2pgsql -s 102722 N2110610 base.contours > N2110610.sql

psql -U postgres -d CM -f N2110610.sql
Then iterate through the appends (note the -a flag)

shp2pgsql -s 102722 N2110615 base.contours -a > N2110615.sql
psql -U postgres -d CM -f N2110615.sql …

we can script this—I have to confess to having done this in excel with worksheet functions since my Bash shell scripting skills help me none in the Windows world, but a friend told me about FOR…IN…DO in batch scripting, so next time… .

So, loading was easy, and all worked well.  Viewing was fine.  What is my gripe?  B-tree Indexing.  What if I wanted to ensure I had only unique geometric records in my database?  I might initiate the table like this instead:

CREATE TABLE “base”.”contours” (gid serial PRIMARY KEY,
“objectid” int8,
“type” int8,
“elevation” smallint,
“shape_len” numeric);
SELECT AddGeometryColumn(‘base’,'contours’,'the_geom’,'102722′,’MULTILINESTRING’,2);
alter table base.contours add constraint unique_geom unique (the_geom);

Adding a unique constraint to the geometry creates an implicit B-tree index, and therein lies the problem– B-tree indices are limited in size to around 2700 bytes.  Our index exceeds that manyfold.  According to the Postgre documentation “Currently, only B-tree indexes can be declared unique.”

This probably is fine for most datasets, but it really limits the use of a uniqueness constraint on geometry columns of large size.  Since this dataset won’t require updates, my work-around (I hope) will be to load the data in the database, perform a select distinct and dump it into another table to clean up any possible duplicates, but it sure would be nice if I could constrain the table in the first place… .

Povray Viewshed with Trees (finally) (cont.)

Ok, your average terrain-only based viewshed (view is from a road to the southeast, viewshed is in cyan):

Note that based on these estimates, we should be able to see most of the valley walls from this little slice of road.  I don’t buy that.  That section of road is pretty closed in with trees.  Let’s add them:

As you may see, just a little cyan in the southeast mostly along the road.  Here’s where you can see the trees occluding our view of the viewshed a bit.  But, it’s a first good stab, even with such boring results… .  I’ll have to get that section of road in here that peeks out over the gorge.  But first, I’ll have to see if I can get tree stem locations from the lidar, so we can trade up for some veracity in our analysis.  Stay tuned.