Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for the ‘POV-Ray’ Category

Gorilla research in Musanze, Rwanda: Hillshades continued

Posted by smathermather on January 30, 2017

I’ve been working on base cartography for the research area in Rwanda. Unlike here in Cleveland, we have some great topography to work with, so we can leverage that for basemaps. But, it’s such a beautiful landscape, I didn’t want to sell these hillshades short by doing a halfway job, so I’ve been diving deep.


First, some legacy. I read three great blog posts on hillshades. One was from ESRI revealing their “Next Generation Hillshade”. Drawing on Swiss Cartographic Traditions, these are some nice looking hillshades using lighting sources from multiple directions (more on this later):

ESRIs Hillshades

Next, we look to Peter Richardsons’s recent post on Mapzen’s blog regarding terrain simplification.

Terrain Generalization example from Mapzen

Terrain Generalization example from Mapzen

I tried (not nearly as hard as I should have) to understand their code, when I saw a link to Daniel Huffman’s blog post from 2011 on terrain generalization: On Generalization Blending for Shaded Relief.

That’s when I saw the equation:

((Generalized DEM * Weight) + (Detailed DEM * (WeightMax – Weight))) / WeightMax

I’ll let you read these posts, rather than rehashing, but here’s what I did toward adding to them. The gist of Daniel and Peter’s approach is to blend together a high resolution and lower resolution version of the DEM based on a weighting factor. Both use a standard deviation filter to determine where to use the high resolution DEM vs resampled version — if the location is much higher or lower than it’s neighbors, it is considered an important feature, and given detail, otherwise the low resolution version is used (actually, I suspect Mapzen’s approach is only highlighting top features based on their diagrams, but I haven’t dived into the code to verify).

The Need

Excuse the colors, we’ll fix those at the end, but this allows us to simplify something that looks like this:

Multicolor hillshade with no simplification

Multicolor hillshade with no simplification

Into something that looks like this:

Multicolor hillshade with simplification

Multicolor hillshade with simplification

See how the hilltops and valleys remain in place and at full detail, but some of the minor facets of the hillsides are simplified? This is our aim.

I developed a pure GDAL approach for the simplification. It is purely command line, has hardcoded file names, etc, but could be done with a python or other API and turned into a proper function. TL:DR: this is not yet refined but quite effective.

Landscape Position

If you’ve been following my blog for a while, you may recall a series of blog posts on determining landscape position using gdal.

Landscape position as calculated on a DEM

Landcape position

This, with small modification, is a perfect tool for determining where to retain DEM features and where to generalize. The one modification is to calculate standard deviation from our simple difference data.

The Tools


Back to those ugly colors on my hillshade version of the map. They go deeper than just color choice — it’s hard not to get a metallic look to digital hillshades. We see it in ESRI’s venerable map and in Mapbox’s Outdoor style. Mapzen may have avoided it by muting the multiple-light approach that ESRI lauds and Mapbox uses — I’m not sure.

HDRI (Shading)

To avoid this with our process (HT Brandon Garmin) I am using HDRI environment mapping for my lighting scheme. This allows for more complicated and realistic lighting that is pleasing to the eye and easy to interpret. Anyone who has followed me for long enough knows where this is going: straight to Pov-Ray… :


The results? Stunning (am I allowed to say that?):

Example of simplified and HDRI rendered hillshade.

Example of simplified and HDRI rendered hillshade.

Hillshade over the hills of Rwanda

The color is very simple here, as we’ll be overlaying data. Please stay tuned.

Posted in GDAL, Gorillas, Karisoke, Optics, POV-Ray, QGIS, R | Tagged: , , , , , , , | Leave a Comment »

Landscape Position using GDAL — PT 3

Posted by smathermather on November 25, 2014

More landscape position pictures — just showing riparianess. See also


valleyz3 valleyz2 valleyz1 valleyz

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | Leave a Comment »

Landscape Position using GDAL — PT 2

Posted by smathermather on November 24, 2014

Just pretty pics today of estimated riparianess. If you prefer a bit of code, see previous post

valleys_improved valleys_improved_1 valleys_improved_2 valleys_improved_3 valleys_improved_3.1 valleys_improved_4 valleys_improved_5

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | 4 Comments »

Landscape Position using GDAL

Posted by smathermather on November 22, 2014

Hat tip again to Seth Fitzsimmons. I’ve been looking for a good, easy to use smoothing algorithm for rasters. Preferably something so easy, I don’t even need to write a little python, and so efficient I can run it on 30GB+ datasets and have it complete before I get distracted again by the next shiny project (a few hours).

Seth’s solution? Downsample to a low resolution using GDAL, then sample back up to a higher resolution in order to smooth the raster. My innovation to his approach? Use Lanczos resampling to keep location static, and get a great smooth model:

Unsmoothed DEM

Unsmoothed DEM

Smoothed DEM

Smoothed DEM

Code to do this in gdal follows. “-tr” sets our resamping resolution, “-r lanczos” sets our resampling algorithm, and the “-co” flags are not strictly necessary, but I’ve got a 30GB dataset, so it helps to chop up the inside of the TIFF in little squares to optimize subsequent processing.

gdalwarp -tr 50 50 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem.tif oh_leap_dem_50.tif
gdalwarp -tr 10 50 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem_50.tif oh_leap_dem_10-50.tif

At first this excited me for cartographic reasons. We can use this to simplify contours, and then use simplified contours at different zoom levels for maps:

But, we can also use this for analyses. For example, if we difference these smoothed images with our original digital elevation model, we get a measurement of local elevation difference, the first step in establishing where valleys, ridges, and other land forms are.

# Resample to lower resolution
gdalwarp -tr 328.0523587211646 328.0523587211646 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem.tif oh_leap_dem_328.tif
# Upsample again to get nicely smoothed data
gdalwarp -tr 3.048293887897243 3.048293887897243 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem_328.tif oh_leap_dem_3-328.tif
# Merge two datasets together into single image as separate bands to ensure they are the same dimensions
# (gdal_calc, as a wrapper for numpy requires this)
gdal_merge -separate -o oh_leap_dem_3-328_m.tif oh_leap_dem.tif oh_leap_dem_3-328.tif
# And now we'll use gdal_calc to difference our elevation model with the smoothed one to get relative elevation 
gdal_calc -A oh_leap_dem_3-328_m.tif -B oh_leap_dem_3-328_m.tif --A_band=1 --B_band=2 --outfile=oh_leap_dem_lp_328.tif --calc="A-B"

So, if we want a good proxy for riparian zones, we can use a technique like this, instead of buffering our streams and rivers a fixed distance (in this case, I’ve used 4 different distances:

Map of landscape position estimated valleys in Cuyahoga County, Ohio

Map of landscape position estimated valleys in Cuyahoga County, Ohio

Pretty snazzy riparian finder. It seems to also find upland headwater wetlands (most of them historic and drained for Cuyahoga County). I am now running on 4 million acres of Ohio at a 10ft (~3 meter) resolution. It’s that efficient.

Addendum: It also finds escarpment edges, like the Portage Escarpment in the above, so it is a mix of a few landforms. Darn handy nonetheless.

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | 2 Comments »

pov-viewshed on GitHub

Posted by smathermather on July 21, 2013

GitHub week for me:

More to come, but in the meantime, enjoy.

Viewshed with Trees
Viewshed, terrain only

Viewshed Summer

Viewshed Winter

Posted in 3D, Analysis, github, Optics, POV-Ray, Viewshed | Tagged: , | 1 Comment »

shape-pov on GitHub

Posted by smathermather on July 20, 2013

I am working on a project, well, I guess just read the README from the repo (a project which will require a creative rename in the future):

ATM, this is a script to write shapefiles to PovRay docs for rendering. Long-term, will be a map tile renderer with 3D support and more.

Whats more? Well, while compositing in TileMill and all the Mapnik derivatives is cool, this tool will handle true material properties, including subsurface scattering, glass like tranparency, diffusion, refraction, and a range of other sophisticated optical effects. The initial focus will be on optical alternatives to compositing rather than 3D.

Stay tuned… .

BTW, will likely be ported from a BASH script to Postgres/PostGIS functions in the future. That and a little middleware glue would provide the interface to the web. This is just the (not yet functional) rapid development version.

Posted in 3D, Analysis, github, Optics, POV-Ray | Leave a Comment »

Tree Interpolator on GitHub

Posted by smathermather on July 20, 2013

I’ve shared bits and pieces of tree interpolation code (3D tree canopies from LiDAR) on this blog for a few years. I finally got things together enough to post some of this code on GitHub. Currently it has a LasTools dependency, so you’ll have to watch for that license gotcha, but I’ll be phasing out that dependency in future versions.

In the future, I’ll do the same with my viewshed work, and some other projects.


Posted in 3D, Analysis, Optics, POV-Ray | Leave a Comment »

Way beyond red-dot fever– bees hives and overlapping home ranges

Posted by smathermather on April 6, 2013

Last week, post Boston Code Sprint, I spent a couple of hours playing with bee data, specifically bee keeper data for Norfolk County Massachusetts.

My friend Eric (a bee keeper of 4 hives in said county) says that worker bees can fly as far as 3 miles for nectar, but after that approximate limit, they hit diminishing returns relative to the calories they burn.

Proximity of bee hives is important for reasons of competition (although bees are quite friendly with each other on the flower) but also disease transfer, which in this age of not-so-easy bee keeping, is an important factor.  So we set out to map bees.

With a list of addresses, QGIS, and the mmqgis extension (see Plugins:Fetch Python Plugins)  which uses Google services for address lookup, we get a nice little set of points for beekeepers in Norfolk County.

01 points

Having a good context map for any such set of points is important.  This is the part of the project I worried most about.  Fortunately, we also have the OpenLayers Plugin, which allows us to add a number of web maps as background (BTW, when we want print quality maps down the road, this will cause us suffering, but for now it is a nice stand-in for a good base layer).  We decided the colorless Stamen Toner map was perfect for our purposes.

02 points toner

Now, let’s buffer these points– give those bees 1.5 miles of free reign, and see where the overlaps in range are (we could make these 3 mile radii, but you get the point):

03 buffer

The map tells us something of overlapping ranges, etc., but it seems like we could use something like we could do something more sophisticated.  Anyone who has followed my blog for a while will not be surprised by the employment of Voronoi polygons… .

04 voronoi

Eric, having suffered an explanation of different ways to calculate Voronoi earlier in the week, remembered the relationship between buffers and Voronoi, and asked if we could combine them.   In this case, picture soap bubbles of your youth.  Blow a soap bubble for each bee point, and the walls of bubbles in between form Voronoi polygons.  In our case, we use cones instead of soap bubbles, but after a couple hours of playing in PovRay, we made some code like this:

#version 3.6;

// Pov Includes
#include “”
#include “”
#include “”

#declare Camera_Size = 500000; //freekin huge

background {color <0, 0, 0>}

//Set the center of image
#declare scene_center_x=759496;
#declare scene_center_y=2861454;

#declare Camera_Location = <scene_center_x,5,scene_center_y> ;
#declare Camera_Lookat = <scene_center_x,0,scene_center_y> ;

// Use orthographic camera for true plan view
camera {
location Camera_Location
look_at Camera_Lookat
right Camera_Size*x
up Camera_Size*y

// Union all the cones together into one object

#declare coner = union {
cone { <0,0,0>,15840,<0,5,0>,0

} };

union {

#declare LastIndex = dimension_size(cone_coords, 1)-2;
#declare Index = 0;
#while(Index <= LastIndex)
object {
translate cone_coords[Index]
#declare Index = Index + 1;

// Pigment cones according to distance from camera
pigment {
gradient x color_map {
[0 color Blue ] // Blue to orange will help us employ blue as the alpha values later
[1 color Orange ]
scale <vlength(Camera_Location-Camera_Lookat),1,1>
Reorient_Trans(x, Camera_Lookat-Camera_Location)
translate Camera_Location
finish {ambient 1}

We get a funky result (yes, this is the result I hoped for):

05 raw pov blue 06 dialog 1

A little transparency, alternate QGIS blending mode, and using the blue band as a transparency band give us a nice output:

07 final 1

We’ll add in those Voronoi boundaries to emphasize the edge of these areas:

07 final

And since this is really only meant for Norfolk County bee keepers, and doesn’t represent folks outside that area, we’ll use a mask with a blending mode of “Subtract” to really offset this map:

08 mask norfolk 09 norfolk masked

I will be interested to see if this helps them look at disease control and spread.  In case you are wondering, Eric is luck enough to have an almost dedicated area for his bees, which includes a huge wetland complex (think: lots of flowers) that is (in his words) all his own.  Well, all his bees own, anyway.

Posted in Analysis, Optics, Other, POV-Ray, QGIS | Tagged: , , , | 2 Comments »

3D Renderers– a woefully incomplete survey

Posted by smathermather on February 4, 2013

Recent exploration into PostGIS 3D tools has reminded me of my love of renderers.  First, there’s PovRay.  PovRay is and always will be my first love.  I’ve done so many fun things with it.

Screen shot 2013-02-04 at 8.10.20 PM

In short, PovRay is a mathematicians renderer– mathematical primitives, advanced ray tracing tools, and a C-like control language make for an abstract but powerful environment.  I like to use it with geographic data, arranged in a planar coordinate system with an orthogonal camera, boom– map.  Check out my povray blog posts to see some of the fun I have had there:

Sadly, the project has been pretty quiet of late (but never goes away, which is nice).  But there are other fun renderer kicking around.  For one, there’s the Physically Based Renderer.  This isn’t just an open source 3D renderer, but a book that teaches you about writing a physically-based renderer, which is really cool.  So far, I haven’t ponied up the $70, nay $61 for the book at Amazon, but one day I will.  What I like about PBRT, it can render arbitrary spectral bands, making it a great remote sensing simulation tool (PovRay can be patched to do this too).

Screen shot 2013-02-04 at 8.03.38 PM

Finally, a renderer that has caught my attention is Mitsuba, which is inspired by PBRT.

Screen shot 2013-02-04 at 8.06.34 PM


Posted in 3D, Other, POV-Ray | Leave a Comment »

Cartographically Pleasing Vegetation Boundaries– with PovRay, PostGIS, and LiDAR

Posted by smathermather on May 18, 2012

Just pretty pictures today.

Posted in Cartography, Database, Optics, Other, PostGIS, PostgreSQL, POV-Ray | Tagged: , , , , | Leave a Comment »