Hillshade (cont.)

In general, I like the results from the hillshade.  One problem that I forgot to take into account– hillshade algorithms typically leave a 1-pixel border along the edge where calculations cannot take place… .  For this Ohio dataset, that means that every 5000 feet or so (the x and y dimension of our DEM chunks), we have a 2-pixel wide line, resulting in a subtle and interesting grid.  Now this is a state plane grid, so I might be able to pretend it’s intentional… .  Maybe I need to mosaic the DEM and then compute hillshade.

Hillshade using PerryGeo’s ported GRASS utility

Using the same LiDAR DEM from which we generated the contours, we can create hillshade tifs using  http://www.perrygeo.net/wordpress/?p=7.  It compiles easily on a mac, probably even easier on a Linux machine following his directions.  Then another simple loop:

#!/bin/bash
x=0
for f in $( ls *.txt); do
 x=`expr $x + 1`
 echo $x $f
 hillshade $f $f.tif
done

I’d like this all in a single tif:

gdal_merge.py *.tif

or in my case with such a large area:

gdal_merge.py -co "BIGTIFF=YES" *.tif

Then overviews:

gdaladdo -r average out.tif 2 4 8 16 32 64 128 256 512 1024

Spreadsheet Fun– Excel

Suppose you had a list of addresses coded like this:

00014 SOUTH ST

That you’d prefer coded like this:

14 SOUTH ST

If you were addicted to modifying things programmatically in Excel, like I quite guiltily do, you could do this in MS Excel.

=(MID(A1,1,FIND(" ",A1)-1))*1 & MID(A1,FIND(" ",A1),LEN(A1))

Everything before the ampersand (&) uses the =mid() function to grab everything up until the first space and multiplies it by one.  Excel automatically does a text to number conversion if you do any mathematical manipulations, so this converts our text to a number, therefore trimming off the spare zeros.  You could use the =value() worksheet function here, but why bother?  After the ampersand, we use the mid function to grab everything after and including the first space.  The concatenation of the two is our cleaned up address.

I know.  Not really free and open source GIS at all… .

GDAL Contours (cont.)

Well, I made some mistakes in the last post, not the least of which is I used the wrong flag for creating an attribute field with elevation.  What follows is a little more sophisticated.  It takes us from a series of DEM tiles from which I want 2-foot and 5-foot contours (using gdal_contour), and then dumps those shapefiles into PostgreSQL using shp2pgsql.

First we prep the database.  Again, I used pre-built binaries from Kyng Chaos.

I created a database called “Vinton” using my favorite PostgreSQL management tool, PGAdminIII (actually, it’s the only one I’ve ever tried…).

So, let’s make Vinton spatially aware:

sudo su - postgres -c '/usr/local/pgsql/bin/createlang plpgsql Vinton'
sudo su - postgres -c '/usr/local/pgsql/bin/psql -d Vinton -f /usr/local/pgsql/share/contrib/postgis.sql'
sudo su - postgres -c '/usr/local/pgsql/bin/psql -d Vinton -f /usr/local/pgsql/share/contrib/spatial_ref_sys.sql'

Vinton county is in Ohio State Plane South.  Our dataset is a foot version of the State Plane South HARN.  Looking it up on SpatialReference.org, I get a description not unlike the PostGIS insert statement below, but I changed the EPSG number from 102323 to 1023236 (arbitrarily), and changed the units to feet (from the default for HARN, which is meters):

UNIT["Foot_US",0.30480060960121924]

so that the full insert is:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 91023236, 'esri', 1023236, '+proj=lcc +lat_1=38.73333333333333 +lat_2=40.03333333333333 +lat_0=38 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs ', 'PROJCS["NAD_1983_HARN_StatePlane_Ohio_South_FIPS_3402",GEOGCS["GCS_North_American_1983_HARN",DATUM["NAD83_High_Accuracy_Regional_Network",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",600000],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-82.5],PARAMETER["Standard_Parallel_1",38.73333333333333],PARAMETER["Standard_Parallel_2",40.03333333333333],PARAMETER["Latitude_Of_Origin",38],UNIT["Foot_US",0.30480060960121924],AUTHORITY["EPSG","1023236"]]');

Now let’s create some tables to house our future contour data:

CREATE TABLE contours_2
(
 gid serial NOT NULL,
 id integer,
 elev numeric,
 the_geom geometry,
 CONSTRAINT contours_2_pkey PRIMARY KEY (gid),
 CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
 CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL),
 CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = 91023236)
)
WITH (OIDS=FALSE);
ALTER TABLE contours_2 OWNER TO postgres;

CREATE TABLE contours_5
(
 gid serial NOT NULL,
 id integer,
 elev numeric,
 the_geom geometry,
 CONSTRAINT contours_5_pkey PRIMARY KEY (gid),
 CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
 CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL),
 CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = 91023236)
)
WITH (OIDS=FALSE);
ALTER TABLE contours_5 OWNER TO postgres;

And finally, we can loop through all of our files, and dump the results first to a shapefile, then pipe that from shp2pgsql directly into the database:

#!/bin/bash
x=0
for f in $( ls *.txt); do
x=`expr $x + 1`
echo $x $f
echo Two Foot:
gdal_contour -i 2 -a elev $f $f.shp
/usr/local/pgsql/bin/shp2pgsql -s 91023236 -a $f.shp contours_2 | /usr/local/pgsql/bin/psql -U postgres -d Vinton

echo Five Foot:
gdal_contour -i 5 -a elev $f 5foot/$f.shp
/usr/local/pgsql/bin/shp2pgsql -s 91023236 -a 5foot/$f.shp contours_5 | /usr/local/pgsql/bin/psql -U postgres -d Vinton
done

Now, I’ll probably convert all those numeric contour fields to integer with a cast statement, but tonight, I let the little laptop do the work… .

GDAL Contours

Just another quick vignette.  From the Ohio Statewide Imagery Program (OSIP) there is a 1-meter DEM for all of Ohio.  To get contours from this dataset, one approach is to use GDAL tools, i.e. gdal_contours.  As I’m working on a Mac today, I used Kyng Chaos pre-compiled Frameworks:

http://www.kyngchaos.com/software:frameworks

Then I needed to update my path variable in the BASH shell:

export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH

Now we can loop through, again, inelegantly (we’ll get unwanted extensions, i.e. *.txt.shp) and create 1-ft contours for a whole county.  Why 1-foot (the OSIP dataset is meant for 2-foot contours).  Well, 1-foot so that later if I want 5-foot contours, I can have them… .

#!/bin/bash
x=0
for f in $( ls *.txt); do
 x=`expr $x + 1`
# echo number of file we are on + name of file
 echo $x $f
# perform analysis, output to *.shp
 gdal_contour -i 1 -nln elev $f $f.shp
done

And the stdout and stderr:

1 S1890435.txt
0...10...20...30...40...50...60...70...80...90...100 - done.
2 S1890440.txt
0...10...20...30...40...50...60...70...80...90...100 - done.
3 S1890445.txt
0...10...20...30...40...50...60...70...80...90...100 - done.
4 S1890450.txt
0...10...20...30...40...50...60..

etc...

LiDAR data in BASH using libLAS

At home, I work on a Mac, so when I want to do work here, my scripting language changes.  In this case, I prefer BASH as my control script.  I downloaded libLAS and did the standard Unix ./configure, make, sudo make install dance in the decompressed libLAS directory

./configure
make
sudo make install

Now I have some tools to manipulate LiDAR las files at home.  To convert all my las files to text:

#!/bin/bash
for f in $( ls *.las); do
 las2txt --parse xyzinrcpM --header pound $f $f.txt
done

There are more elegant ways to implement this, but it gets the job done… .  Now I have a few hundred text files with all the relevant LiDAR info in them for manipulation.  Since working with the whole enchilada in PostGIS failed, I think I’ll approach it by loading one section at a time, analyzing, and removing those data from the table.  More to come… .

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/