Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for January, 2010

Mapping places unknown– free global datasets and FOSS GIS are a great combo

Posted by smathermather on January 31, 2010

I wanted to put together a quick and dirty map of a biological reserve in Ecuador, sort of a laptop exploration of a place quite distant.  At first, I thought I’d use Shuttle Radar Topography Misssion data to get the elevation information.  Then I discovered the ASTER Global DEM which is 30m resolution for the whole world.  Wow.  Cool cool data. (I used the Japanese site to download, as I had trouble getting the US site to load).

So, first I downloaded and extracted the zipped geotiffs that come from the ASTER Global DEM.  I extracted the DEMs into the same directory, and mosaicked them into a single DEM:

gdalwarp AST*.tif ecuador_aster_dem.tif

Then I created a hillshade version of the DEM for mapping:

gdaldem hillshade ecuador_aster_dem.tif ecuador_aster_hillshade.tif

and then 100 and 50 meter contours for the site:

gdal_contour -a elev -3d -i 100 ecuador_aster_dem.tif ecuador_aster_contour_100.shp
gdal_contour -a elev -3d -i 50 ecuador_aster_dem.tif ecuador_aster_contour_50.shp

I loaded it all into Quantum GIS, symbolized it, and added some OpenStreetMap data.  To do this, I navigated to the location in OSM, chose the download tab, and chose OpenStreetMap XML Data as my output.  Using the OSM plug-in in QGIS, I imported the file, and walla– a nearly functional map.

Posted in GDAL | Tagged: , , , | Leave a Comment »

Clipping Data w/ PostGIS continued

Posted by smathermather on January 29, 2010

It finished, yay!  A little Friday 5PM exuberance, I’m afraid.  See previous post for explanation.

Here is the clipped version:

Posted in Database, PostGIS, SQL | Tagged: , , , , , | Leave a Comment »

Clipping Data

Posted by smathermather on January 28, 2010

I loaded a set of contours into the database that are really nice, up-to-date LiDAR and breakline derived contours.  The are engineering grade aerial contours, but they are a very big and complex dataset.  So I’ve done what I had hoped was the hard part, using shp2pgsql, I’ve converted them to PostGIS insert statements and dumped them into the database only to find that some sections are corrupt.  I’m hoping to use the function ST_Difference to clip out the bad areas, so I can insert them in again.

There are two problems in the contour dataset– some areas are duplicated records, and some areas are missing or otherwise corrupted areas.  I’ve got a polygon file (in purple below) that delineates where the problem areas are, so I hope to use the function ST_Difference to remove these areas so I can repopulate with good data.

CREATE TABLE public.cuy_contours_2_clip
 AS
 SELECT elevation,
       ST_Difference(a.the_geom, b.the_geom) AS the_geom
    FROM base.cuy_contours_2_1 as a, 
         public.contour_clip_bounds as b;

Let’s see if PostGIS can handle this big of a dataset… .

Posted in Database, PostGIS, SQL | Tagged: , , , , , | 2 Comments »

Viewshed Analyses in Povray– Final Image

Posted by smathermather on January 28, 2010

I completed this project a long time ago and have not blogged on it in over a year.  None-the-less, I realized that I hadn’t posted any final images for the viewshed analysis with Povray.

So here is the summer viewshed analysis.  This is rendered with about (I think) 150,000 trees, each with 100,000+ leaves.  The cell tower here is in red.  The landscape colored in cyan where the tower can be seen from.  A property boundary shows up in green:

And here is the same viewshed but a winter scene.  Since this is a very flat wetland landscape, it’s important to understand the visual penetration of the cell tower with the vegetation, with and without leaves:

Posted in POV-Ray, Viewshed | Tagged: , , , , | 1 Comment »

IMCORR– using image correlation to georeference an image

Posted by smathermather on January 23, 2010

IMCORR is a package distributed by the National Snow and Ice Data Center (did you know we have one of those?), or NSIDC, that performs image cross correlation between two images using a comparison between a moving image chip in each image.  It measures the displacement (in pixels) between the objects found in the two images, and writes that out to a text file.  It is used to calculate how quickly glaciers are moving by calculating displacements.  If you have an image of a glacier from two different time periods, this little program will calculate how far crevasse zones and other discernible features have moved between two images. I know I know.  I have been talking about trees, and now ice.  Bear with me.

IMCORR can also be used for simple image georeferencing, say when we don’t want to collect our own tie points.  We can look for correlations between two images, and tabulate those as points, use them with a thin plate spline, in say gdalwarp, and thus correct for distortions in one dataset with another one.

Now, this is all in theory, and I haven’t tested it yet.  But I’m thinking this may be a really useful technique, especially because we already have well referenced high-resolution aerials for so many areas.  I have an effective 3″ color infrared (CIR) dataset for an entire county kicking around somewhere that has been referenced to frame on center, but with no terrain correction or other optical distortion corrections applied.  I’m hoping to cross correlate these with a 6″ georeferenced aerial in order to make these data useable.  3″ CIR could be very useful for all sorts of fun applications… .

Posted in Other | Tagged: , , , , , , | 1 Comment »

Digital Surface Model– a whole forest and more part Deux

Posted by smathermather on January 18, 2010

Just another shot of our digital surface model of the forest.  This time, I rendered it with all of the trees at base elevation 0 (zero), placed the orthographic camera at height 150, so our 16-bit range we’re rendering to should scale 0-150 feet to 0-65535 values, giving us a precision of 0.002 feet.  The original LiDAR data has a precision of 0.01 feet, so we’re not throwing data away with the integerization.  By the way, we don’t have enough LiDAR points to really derive canopy shape– we’ve invented the detailed canopy shape by using a tree shaped interpolation technique (if we can quite liberally call it interpolation– I’d hate to see the actual function written out).

Anyway, the pretty pic– the darker the green, the taller the tree– we max out around 120 feet in most areas.  We do find at least one 149 foot tall tree, but I’m guessing that’s on a steep slope, and is an over estimate due to a tree leaning of spreading over a cliff, or some such scenario.  Of course, if it does exist, I should go look for it.  I did measure one in the field last year that exceeded 130 feet, and I couldn’t spot the top, so it may have been much taller.

TreeDEM

Oh, and Mac bragging rights here– I couldn’t render the full 30,000 x 30,000 pixel scene on Windows XP (although Vista may have handled it), but it rendered just fine in 35 minutes on my Mac laptop.

Posted in POV-Ray | Tagged: , , , , , | Leave a Comment »

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.

Posted in POV-Ray | Tagged: , , , , | Leave a Comment »

Topographic Position Index and Ecological Land Type (warning completely unrefined not quite Geologic dribble– with bad maps :) …)

Posted by smathermather on January 10, 2010

Warning.  What follows is somewhat informed, but I’m no geologist.  I just play one on wordpress.

Understanding the basic underlying geology and associated topography plus site history helps us achieve a basic understanding of a sites ecological potential.  At the most basic level, we expect different wildlife and vegetation dynamics in a floodplain vs. a mountain ridge.  Classification of digital elevation models can be done with a concept called topographic position index (of which Jensen Enterprises has a pretty good explanation).

I’ve been playing with TPI for two areas in the Allegheny Plateau in Ohio– one in an unglaciated portion, and one in a glaciated portion.  The Allegheny Plateau developed as the once very tall Appalachians wore down to the nubs they are today.  That plateau has since eroded into little hill-like remains, with stream valleys redefining a hill and hollow landscape, and most of the ridge-tops being near the same elevation.  In the glaciated portion, this is overlayed with glacial till, a re-flattening of these “hilltops”, and the definition of gorges in places where sub-glacial rivers subdivided the terrain.

The TPI attempt I did for the unglaciated portion roughly defines the ridge tops in red, stream valleys and slope bases in blue, and slopes as yellow/yellow green.  The map that follows is the overlay of TPI performed in gdaldem at multiple scales, but we are somewhat limited in what we can do with gdaldem, because we can’t control our window size.

Then I performed a proper TPI analysis at 80 feet and 300 feet for a glaciated portion of the plateau in Northern Ohio, following Jensen’s guidlines.  What follows is a photo of a printed map, cause I forgot to bring an jpeg of pdf home… .  Here, upland flat portions are dark green.  They are adjacent to light blue areas that are also upland, but exceed 5% slopes.  We can see upland areas adjacent to cliff areas as well as peaks in red, lower areas of high relief in orange, gorge valleys in dark blue, U-shaped valleys in a lighter blue (stream polygons and lines are also overlayed, adding to color confusion– hence the warning in the title).

What intrigues me about this map are the areas around the streams just northwest of the center of the map.  Capturing the gorges and upland areas successfully was not a surprise.  What is surprising is that we capture and map the areas that seem to me likely contribute directly to surface flow into the streams as distinct areas from the upland plains.  Next, I think I’ll compare this analysis with the soil map for the area, as well as the detailed vegetation cover we have for the area.

Posted in Ecology, Landscape Position, Optics | Tagged: , , , , , , , | 3 Comments »

Modeling (relative) Sub-Canopy Biophysical Variables with PovRay

Posted by smathermather on January 9, 2010

One of my goals in modeling canopy height is to estimate sub-canopy biophysical variables that give indications of habitat suitability.  For example, suppose we have a rare plant that is canopy gap dependent.  How do we model it’s distribution?  With a statewide LiDAR dataset from which we can derive canopy location and height, we can identify canopy gaps.  While I haven’t gotten to the modeling phase yet, my maps of canopy boundaries have been really neat to field verify– it is possible for me to identify where I am in the forest by looking at at the size and shape of the canopy gap I am near, relative to other features, such as streams and topography.

Another use of this facility is to model urban farming opportunities– sunlight (site aspect and shading) and water being the two primary driving factors for choosing urban agriculture locations (soil and/or soil ammendments are typically hauled in if the soil is inadequate).

PovRay isn’t the best tool for doing this analysis (we can only model relative shadowing, not absolute physical solar incidence), but a Swiss Army knife isn’t the best bottle opener either.  I’ve got one of each in my pocket, so I use them.  PovRay does have a sun position model controlled by latitude, year, time of year, day, hour, etc.

Shadow length at 12:34 every day for a year:

#include "sunpos.inc"

global_settings {assumed_gamma 1}

camera {
 location  <3.0, 5.0, -5.0>
 look_at   <0.40, 1.8, 0.0>
 angle 40
}

#declare Year= 1998;
#declare Month= 1;
#declare Day= 1;
#declare Hour= 12;
#declare Minute= 34;
#declare Lstm= 15;
#declare LONG= 6.9;
#declare LAT=42;

#declare LastIndex = 365;

#declare Index = 1;
#while(Index <= LastIndex)

 //Put in the suns
 light_source {
 SunPos(Year, Month, Day, Hour*LastIndex, Minute, Lstm, LAT, LONG)
 rgb 1/LastIndex
 }

#declare Index = Index + 1;
#end

background {rgb 0}
#include "Linden.inc"
 object { TREE
 scale 4
 }
plane {
 <0, 20.0, 0.0>, 0.0
 material {
 texture {
 pigment {
 rgb <0.700023, 0.700023, 0.700023>
 }
 normal {
 brick 0.1 //amount
 }
 finish {
 diffuse 0.6
 brilliance 1.0
 }
 }
 interior {
 ior 1.3
 }
 }
}

12:34, January 1st:

And the same tree and shadow, 179 days later:

And the summation of all shadows from above:

#include "sunpos.inc"

global_settings {assumed_gamma 1}

camera {
 orthographic
 location  <0, 5.0, -2>
 look_at   <0.40, 1.8, 0.0>
 angle 80
}

#declare Year= 1998;
#declare Month= 1;
#declare Day= 1;
#declare Hour= 12;
#declare Minute= 34;
#declare Lstm= 15;
#declare LONG= 6.9;
#declare LAT=42;

#declare LastIndex = 365;

#declare Index = 1;
#while(Index <= LastIndex)

 //Put in the suns
 light_source {
 SunPos(Year, Month, Day, Hour*LastIndex, Minute, Lstm, LAT, LONG)
 rgb 1/LastIndex
 }

#declare Index = Index + 1;
#end

background {rgb 0}

#include "Linden.inc"
 object { TREE
 scale 4
 no_image
 }

plane {
 <0, 20.0, 0.0>, 0.0
 material {
 texture {
 pigment {
 rgb <.7,.7,.7>
 }
 normal {
 brick 0.1 //amount
 }
 finish {
 diffuse 0.6
 brilliance 1.0
 }
 }
 interior {
 ior 1.3
 }
 }
}

Posted in Ecology, Optics, POV-Ray | Tagged: , , , | 1 Comment »

Optimizing Pov-Ray include files– smaller trees with the same effect

Posted by smathermather on January 8, 2010

One limitation I’ve run into in rendering large scenes of trees in PovRay, whether for fun for work is memory usage.  I’ve found that about a half million trees is about as many as I can place on a 32-bit system before running out of memory.  Sometime this year I’ll switch to 64-bits, but in the mean time, I hacked together a solution for reducing my tree include file sizes to save memory.

My include file is called tree.inc.  It has the following format (generated by PovTree):

#declare FOLIAGE = mesh {
triangle{<0.18162394, -0.24128051, 0.030293489>, <0.18287745, -0.23965627, 0.026677>, <0.18274091, -0.23925833, 0.033750653>}
triangle{<0.18287745, -0.23965627, 0.026677>, <0.18274091, -0.23925833, 0.033750653>, <0.18468907, -0.23614006, 0.0350326>}
...
...

#declare TREE = union {
object{FOLIAGE}
}

You’ll notice 8 decimal places for those coordinates.  We probably don’t need that level of precision.  The code that follows is messy, but it does the trick– using printf within awk to reduce the total number of decimal points (3 decimal places in this case).  If memory serves me, we strip all the extra brackets, carrots, and commas with sed, format with awk as a 3 decimal place number, put our brackets, carrots, and commas back where they belong, and the trim out any extra characters with sed before piping to a new include file, tree1.inc.  Oh, and I think I just manually added back on the FOLIAGE declaration, etc..

more tree.inc | grep "triangle" | \
sed -e 's/{/ /g' -e 's/</ /g' -e 's/>/ /g' -e 's/}/ /g' -e 's/,/ /g'| \
awk '{ printf "%0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f \n", $2, $3, $4, $5, $6, $7, $8, $9, $10};' | \
awk '{ print "triangle {<", $1, ",", $2, ",", $3, ">, <", $4, ",", $5, ",", $6, ">, <", $7, ",", $8, ",", $9, ">}"    };' | \
sed -e 's/ //g' -e 's/0\./\./g'  > tree1.inc

resulting in:

#declare FOLIAGE = mesh {
triangle{<.182,-.241,.030>,<.183,-.240,.027>,<.183,-.239,.034>}
triangle{<.183,-.240,.027>,<.183,-.239,.034>,<.185,-.236,.035>}
...
...

#declare TREE = union {
object{FOLIAGE}
}

Tree image before:

Tree image after:

What’s the difference in size for the include file?

-rw-r--r--  1 user  group    13M Nov  6 21:52 tree.inc
-rw-r--r--  1 user  group   5.6M Nov  8 00:40 tree1.inc

How far can we go with this?  Well, here’s two decimal places for our coordinates:

It’s starting to get a little chunky here, but it might work for a lot of purposes.  However, we’re reaching the limit of our optimization:

-rw-r--r--  1 user  group   5.2M Nov  8 00:48 tree4.inc

Posted in POV-Ray | Tagged: , | Leave a Comment »