Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Digital Surface Model– a whole forest and then some

Posted by smathermather on January 12, 2010

Let’s start putting some of the pieces together.  Earlier, I presented a method for deriving a digital surface model from LiDAR data of a forest canopy, using a tree-shaped interpolator that scales the individual tree canopy shapes according to the height of the LiDAR return from the ground.  Now I present the whole project, forest and all.   I wanted to render my DSM at a foot resolution, but that meant rendering a 30,000 x 30,000 pixel image.  PovRay would have none of that on my 32-bit system, especially since it was also rendering about a half million trees that it had to keep in memory.  So, I compromised, and tiled out the rendering process by rendering 17 5,000 x 5,000 images for mosaicking back together (not all 36 tiles had data in them).  My include files include coordinates for tree locations, the tree height include file, and a camera and camera look at include file.

#include "transforms.inc"
#version 3.6;

#include "..//Tree1.inc"
#include "sc_tree_coords.inc"
#include "sc_tree_height.inc"
#include "sc_camera_coords.inc"
#include "sc_lookat_coords.inc"   

#declare Camera_Location = camera_coords[frame_number];
#declare Camera_Lookat   = lookat_coords[frame_number];
#declare Camera_Size = 5000;

camera {
 orthographic
 location Camera_Location
 look_at Camera_Lookat
 right Camera_Size*x
 up Camera_Size*y
}   

background {color <0, 0, 0>}

union {

#declare Rnd_1 = seed (1153);      

#declare LastIndex = dimension_size(tree_coords, 1)-2;
#declare Index = 0;
#while(Index <= LastIndex)

 object  {
 TREE
 scale tree_height[Index]
 rotate <0,rand(Rnd_1)*360,0>
 translate tree_coords[Index]
 }

#declare Index = Index + 1;
#end

 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}
}

sc_tree_coords.inc:

#declare tree_coords = array[435496]
{<2262851.09,1036.09,635014.03> ,
<2262753.22,1036.22,635015.56> ,
<2262850.37,1037.06,635024.82> ,
<2262819.89,1036.61,635024.85> ,...
<2262908.31,1035.49,635024.05> ,
<2263031.45,1036.12,635007.38> ,
<2263020.26,1036.77,635000.08> ,
<2262896.1,1035.72,635014.09> }

sc_tree_height.inc:

#declare tree_height = array[435496]
{18.67,
25.83,
56.87,
41.26,...
55.18,
53.99,
36.2,
64.46}

sc_camera_coords.inc:

#declare camera_coords = array[17]
{<2272500,1500,637500> ,
<2272500,1500,642500> ,
<2272500,1500,647500> ,
<2257500,1500,637500> ,
<2257500,1500,642500> ,
<2247500,1500,627500> ,
<2247500,1500,632500> ,
<2247500,1500,637500> ,
<2252500,1500,637500> ,
<2262500,1500,637500> ,
<2262500,1500,642500> ,
<2262500,1500,647500> ,
<2267500,1500,632500> ,
<2267500,1500,637500> ,
<2267500,1500,642500> ,
<2267500,1500,647500> ,
<2267500,1500,652500> }

sc_lookat_coords.inc:

#declare lookat_coords = array[17]
{<2272500,500,637500> ,
<2272500,500,642500> ,
<2272500,500,647500> ,
<2257500,500,637500> ,
<2257500,500,642500> ,
<2247500,500,627500> ,
<2247500,500,632500> ,
<2247500,500,637500> ,
<2252500,500,637500> ,
<2262500,500,637500> ,
<2262500,500,642500> ,
<2262500,500,647500> ,
<2267500,500,632500> ,
<2267500,500,637500> ,
<2267500,500,642500> ,
<2267500,500,647500> ,
<2267500,500,652500> }

And a very small subset of the output:


This outputs as a 16-bit PNG.  I placed my trees at real heights relative to the DEM I have, so the range of values is absolute height, not height relative to the DEM.  To get height relative to the DEM, I have to subtract the DEM from my output here.

For my camera positions, the difference between the camera location and lookat point is 1000 feet, as my lookat locations were 500 feet in the Z and camera location 1500 in the Z (the full range of elevations in Ohio is about 550 feet to 1550 feet, and the part of Ohio I’m in has neither the highest or lowest elevations). That 1000 feet is scaled by PovRay to 16-bit integer, or 65535.  So to convert the output to real world units, I divide by 65535, multiply by 1000, and then add 500.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: