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

Just pretty pictures today.

Landscape Position: Conclusion? (part 2)

From earlier post:

“I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.”

I didn’t use PostGIS Raster, but the incredibly useful gdal_calc.py to finish out my landscape position calculations.  I will endeavor to post my code to github soon, but the parts and pieces include using gdal and PovRay.  PovRay helps with the sampling (really!) of nearby neighbors in the raster DEM, and gdal does the averaging and differencing of those to get relative landscape position.  I spent some time yesterday scratching my head over how to show all the new landscape position information on a readable and useable map, and after discussion with collegues, decided to use it to divide the world into two categories– riparian & lowland + headwaters & upland (not reflected yet in the labels).  To find out more about landscape position, follow this link, or better yet this one.  (BTW, the green is park land, blue is riparian/lowland/stream networks, purple is the basin boundary).

Riparian map based on landscape position, calculated using PovRay and GDAL.

Landscape Position: Conclusion?

I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.

Dialog– What qualifies as a benchmark? Part 1

Normally, my blog is a bit of a monologue.  It’s not a bad thing, but can be a little lonely.  Every now and then I get a great (and often doubtful) comments, which enhances things considerably.

What follows is some dialog about the LiDAR shootout series, largely between Etienne and Pierre, posted with their permission:

Pierre:

“Etienne, Stephen,

“I really appreciate the benchmark work you are doing to compare different technologies to solve your lidar data processing problem.

“However, if you want your tests to have as much credibility as possible and mostly if you wishes to publish the results on the web I would recommend that you compare horses with horses, dogs with dogs, general vector technologies with general vector technologies, raster technologies with raster technologies, specialized lidar technologies with specialized lidar technologies. Otherwise the dog will always lose against the horses. You have to be clear on that and well explain the differences and the pro’s and con’s between all the tech you will be using in your benchmark. Are you doing the test in ArcGIS using ESRI lidar technology (or general vector one)? Same question for R and PostGIS Raster (in the last case it is clear that the answer is no). In your graph you are comparing very different kind of animals without giving much chances to dogs.

“For every technologies you will also have to test different approaches in order to find the fastest one for this technology in order to compare it with the fastest one in the other technologies. Comparing a slow approach in one tech with the fastest in the other would be a bit unfair. This is what we did for PostGIS Raster. Did you do the same with R and ArcGIS?

“Basically I think it is just unfair to compare ArcGIS, R and PostGIS Raster with LAStools for working with Lidar data. Mostly if you are not using ArcGIS and R packages for dealing with lidar. Teach us you did not use them if you did not. We will all learn something.”

Etienne:

“As you said, the point is to compare different strategies to get the job done. That’s it. It’s not a question of being fair or not, we report facts. Maybe the bottom line is that pgraster shouldn’t be used to extract lidar height. However it can be used for other things. Well, we all learned something.

“It also emphasize the fact that there is no point_cloud type in postgis. I guess this could improve a lot the postgis results as it seems handling point cloud might be really fast. I learned that.

“Each solution has its drawbacks, for ArcGIS I think we’ll all agree about the cost of the licence and furthermore the cost of the lidar plugin. For lastools, it require a licensing for commercial use (which is really low btw compared to ArcGIS + Spatial analyst + LP360). PG raster require a postgres installation and some SQL knowledge, etc. R require to learn a new language, etc.

“Now I must acknowledge, the tests we did were on only 500k points. That could be considered a lot, and not. When dealing with lidar, one could easily expect to get more than a thousand times this quantity. But it was the limit imposed by the LP360 version I used and also a reasonable size for benchmarking. Note that pg raster was really stable in the timing (whatever the size).

“Finally, please note that Martin’s approach is new and he just released the tool a week ago. It’s to the better of my knowledge, the first one the do it.

“Pierre, explain your point better, I’m not sure I get it (or I disagree, as you can see).”

Pierre:
“My point is just that you have to be clear in your post about the fact that LAStool is a technology developed explicitly for dealing with lidar data. PostGIS raster is not and I guess you are not using ESRI lidar technology and I have no idea how you do it with R. The point is to compare potatoes with potatoes.

“If I would be from the ESRI side I would quickly post a comment in the post saying “Did you read this article?” http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Estimating_forest_canopy_density_and_height/00q8000000v3000000/

“A don’t know how you did it in ArcGIS but it seems that they are storing lidar data as multipoints.

“A good benchmark should precise all this.“

Etienne (Pierre > quoted):

> If I would be from the ESRI side I would quickly post a comment in the post saying “Did you read this article?”

> http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Estimating_forest_canopy_density_and_height/00q8000000v3000000/‘
“I would reply: apart from the fact that it uses lidar data, it has nothing to do with the object of the present article But I’d be happy to discover a better way than what I did (for all the tools btw).

“Maybe you should (re-)read the post, it’s written. For R as well. Maybe Steve, you should precise (this is my fault) that for ArcGIS, I’ve used LP360 (evaluation) to load lidar data, which was limited to 500k points. But it wasn’t involved in the computation as I used (as written in the post) «ExtractValuesToPoints».”

Pierre:

> So it’s a good benchmark ?

“I just think it would be more fair to create three test categories:

“1) Lidar technologies (apparently your ArcGIS test should be in this category)

“2) Vector/raster technologies (why not use non lidar techs if lidar techs are so much better?)

“3) Nearest neighbor technique using vector only technology

“I won’t comment further…

Etienne:

“I personally would not encourage these categories as there is some confusion (at least in my mind) on the definition of what is a «LiDAR thechnology». However I think it would be a good idea to nuance the conclusions, which I think will be done in further posts (or not). It’s also pretty clear from the RMSE measures what is the main drawback of each solution (with speed).

“I’m still open for discussion…

Stephen (me):

“Hi Pierre,
     “I agree with you in principle, but the ArcGIS aside, there are some very real opportunities for optimization on point in raster analyses with tiling, indexing, etc.  In other words, I don’t see any practical reasons why PostGIS Raster couldn’t be (nearly) as fast as an nn technique, so long as the tile size is optimized for density of the analysis in question.  The only thing that makes it a dog is its age as a technology.

(moments later after discovering the lengthy dialog…)

“It looks like I missed all the discussion here…

“The reason such benchmarks are interesting is (in part) that they are not fair.  PostGIS isn’t optimized for LiDAR analyses.  PostGIS Raster isn’t optimized for LiDAR analyses.  PostGIS Raster is very young technology which really isn’t fair.  The PostGIS development folks know this, and were interested in seeing my initial benchmarks with nn analyses for LiDAR height (last year) because they are aware of the optimizations that are not yet there for handling point clouds.  Vector and raster analyses (as we’ve run them here) should not in principle or practice substantially different in speed, so long as the query analyzer maximizes the efficiency of the techniques.  Add in multipoint optimizations and other voodoo, and well, that’s a different story.  Another aspect to this is that comparing lastools to the rest is unfair not just because it’s optimized for lidar but because it’s working at a much lower level (file level) than the abstractions a database provides.  The trade-off, off course, is that a database can provide transaction support, standardized language, automated query analyzers, etc. etc. that don’t apply for file level processing.

“What I see here, categories aside, is a need for the other 2-3 parts of this 3-4 part series– the breakdown of which techniques are used, how they’re implemented, and the rationale behind them.  Don’t worry, we’ll get there.  That’s the danger of starting with the punchline, but we’ll build the nuance in.“

 

So, there we are.  Looks like I’m committed to building this out with a little more nuance.    Stay tuned.

LiDAR Shootout! — New Chart, Final Results

Comparison of lidar height calculations

In reviewing the final numbers and charts from Etienne and Pierre, above are the results we see.  The only revision is a moderate increase in speed for the PG Raster query.

Final results in speed for lastools– ~350,000 points per second.  In other words– off-the-charts fast.  And the initial RMSE of ~25 feet was a mistake– it is probably closer to 0.2 feet.

Stay tuned for detailed reviews of these techniques (with code).

 

LiDAR Shootout!

For a couple of months now I’ve been corresponding with Etienne Racine and Pierre Racine out of Montreal Laval University in Quebec City.  They decided to take on the problem of finding the speed and accuracy of a number of different techniques for extracting canopy height from LiDAR data.  They have been kind enough to allow me to post the results here.  This will be a multipart series.  We’ll start with the punchline, and then build out the code behind each of the queries/calculations in turn (category will be “LiDAR Shootout“) with separate posts.

In the hopper for testing were essentially 4 different ways to extract canopy height from a LiDAR point cloud and (in three of the cases) a DEM extracted from the LiDAR point cloud.  The 4 techniques were as follows:

• Use ArcGIS ExtractValuesToPoints.
• Use R using the raster and rgdal packages with raster::extract() function with in-memory results.
• Use postgis (2.0 svn) Raster datatype (formerly WKT Raster)
• And use a simple nearest neighbor approach with ST_DWithin within PostGIS.

The data are from the Ohio Statewide Imagery Program (OSIP) program, run by Ohio Geographically Referenced Information Program (OGRIP).  Ohio is blessed with an excellent 2006-2008 LiDAR dataset and statewide DEM derived from the (effectively) 3.5 foot horizontal posting data (specs may say 5-ft, I can’t remember…).

Speeeed:

So, first to speed, initial results.  Who wins in the speed category?  Measured as points per second (on consistently the same hardware), nearest neighbor wins by a landslide (bigger is better here):

Application:  Points per Second

ArcGIS:  4504
R:  3228
PostGIS Raster:  1451 (maybe as high as 4500)
Postgis nn:  18208

Now, later on, Pierre discovered changes to indexing may help an the EXPLAIN query analyzer optimization which tripled the PostGIS Raster query speed, making it about the same speed as ArcGIS.  More on that later.

Figure removed– to be replaced in another post

Accuracy:

A fast calculation is always better, unless you trade off accuracy in some detrimental way.  With PostGIS NN, we’re not interpolating our LiDAR ground point cloud before calculating the difference, so relative to the other three techniques, we could be introducing some bias/artifacts here, a point Jeff Evans makes here.  Overall error relative to the interpolated solution introduced by using an NN technique on this dataset: 0.268 feet.  Whew!  Sigh of relief for me!  We may spend some time looking at the geographic distribution of that error to see if it’s random or otherwise, but that’s very near the error level for the input data.

Addendum:

8 days ago, Martin Isenburg’s let the lasheight tool drop. Lastools is impressively fast.  Preliminary tests by Etienne: 140450 points per second.

Oy!  A factor of 8 faster than PostGIS NN  And it uses an on-the-fly calculated DEM using delaunay triangulation.  Preliminary results indicate a very different DEM than the one calculated by the state:

RMSE: 25.4227446600656

More research will need done to find out the source of the differences: they are not random.

In other news:

PostGIS will get a true Knn technique in the near future.  Such a project is funded, so stay tuned for faster results:

http://www.postgis.org/pipermail/postgis-devel/2011-June/013931.html

Also, index changes to PostGIS could come down the pike, if funded, that will result in bounding box queries that are 6-times faster, ala:

http://blog.opengeo.org/2011/05/27/pgcon-notes-3/

Between the two optimizations, we could give Lastools a run for its money on speed  .  In the meantime, preliminary RMSE aside, Lastools (the 5th tool not in this test) may win hands down.

Multi-Ring Buffers in PostGIS– Prep for Landuse Analysis

If we want to understand the factors affecting the quality of an ecological resource, adjacency is important. Adjacent land use and the fossil fuel inputs relative to keeping a land use in its current use seem to have strong correlation with, e.g. wetland quality.  I’ll have a deeper post on this with some citations, but this relatively simple process of analysis is called Landscape Development Index (LDI).  For now, just some PostGIS code with pretty little wetland buffer pictures… .  More complete code forthcoming… .

CREATE TABLE public.wetland_bufftest4
	AS

	SELECT wet.gid AS gid_wet, ST_Difference(
		ST_Union(ST_Buffer(wet.the_geom, 100)),
		ST_Union(ST_Buffer(wet.the_geom, 0))
		)
	AS the_geom
	FROM public.wetland wet
	GROUP BY wet.gid

	UNION ALL

	SELECT wet.gid AS gid_wet, ST_Difference(
		ST_Union(ST_Buffer(wet.the_geom, 250)),
		ST_Union(ST_Buffer(wet.the_geom, 100))
		)
	AS the_geom
	FROM public.wetland wet
	GROUP BY wet.gid

	UNION ALL

	SELECT wet.gid AS gid_wet, ST_Difference(
		ST_Union(ST_Buffer(wet.the_geom, 500)),
		ST_Union(ST_Buffer(wet.the_geom, 250))
		)
	AS the_geom
	FROM public.wetland wet
	GROUP BY wet.gid

	UNION ALL

	SELECT wet.gid AS gid_wet, ST_Difference(
		ST_Union(ST_Buffer(wet.the_geom, 1000)),
		ST_Union(ST_Buffer(wet.the_geom, 500))
		)
	AS the_geom
	FROM public.wetland wet
	GROUP BY wet.gid
;

ALTER TABLE public.wetland_bufftest4 ADD COLUMN gid serial;
ALTER TABLE public.wetland_bufftest4 ADD PRIMARY KEY (gid);

Landscape Position and McNab Indices (cont.)

I typed that last one too quickly– too many typos, but my wife says I’m not supposed to revise blogs, but move on… .

So, for clarity, let’s talk a little more about McNab indices.  Field-derived McNab indices are a measure of average angle from the observer to the horizon (mesoscale landform index), or from the observer to another field person a set distance away, e.g. 10m or 18m, or whatever the plot size or settled upon (minor landform index).  Typically these are done in the field at a discreet number of directions, e.g. 4 cardinal directions, or 4 cardinal plus 4 ordinal (NE,NW SE, SW, SE) directions. The landform image in this post and the last is a calculation of mesoscale landform, which is harder to do in a classic GIS than minor landform (I’ll have a follow-up post on minor landform, probably using ArcGIS).

To calculate the values for mesoscale landform computationally, we require calculating the angle to the horizon for a certain number of directions for each point in the landscape.  This has the potential to be computationally intensive and complicated.  If, however, we conceive of the problem differently, as a 3D calculation we can perform in PovRay, arriving at the answer is simplified by the coding already done by the PovRay programmers.

Essentially, McNab mesoscale indices are a field proxy for the steradian of a site.

What does that mean?  Well, essentially, a map of steradians is a proxy for the shading of a white uneven surface on a cloudy day– or how much diffuse light is available to a given site.  Used in combination with site aspect, this is enough information to determine most of of the light conditions of a site, which is why McNab indices in combination with other factors are a good predictor of site productivity, and correlated with different plant communities across the landscape.

How does an uneven white surface shade on a cloudy day? Steeper areas with more of the sky occluded are darker while wide open spaces, like the bottom of river valleys and the tops of ridges and plateaus are brighter.  If you want to witness this effect, look to snow on the ground on a cloudy day (and what a great winter to do it).  (The only difference with snow is subsurface scattering which evens out the light quite a bit.)  You’ll notice the divots in the snow to be darker than the peaks, and the edges of the divots to be darkest of all.

The question then, is how to we compute this within PovRay?  We could use radiance as a global illumination model, but the calculation of inter-reflectivity that is at the core of radiance, while an important real factor in the landscape, would fail to replicate the original mesoscale landform index.  Instead, we set up a series of lights in a dome to illuminate the whole sky sphere, blur them a bit, and call it a day, a technique developed for HDRI rendering by Ben Weston, whose code I borrowed heavily.  The more lights the element of the landscape is exposed to the brighter, and vice versa, essentially making the brightness of an element proportional to the exposure to the sky sphere.  Unlike Ben, I used a simple white sphere to get even lighting.

Rendered at 16-bit resolution, we have a possible range of values from 0-65535.  Let’s assume that a linear transformation of these values will result in values representing the proportion of the sky sphere.  From there, transformation to steradians and then to solid angles in degrees is trivial.  Once it’s solid angles in degrees, it represents the same kind of value as a mesoscale landform index would give us (more later…).

Cost of nearest neighbor search depending on distance

A quick review of costs to search with increasing distances.  Reference this original post for the code being run.

SELECT DISTINCT ON(g1.gid)  g1.gid as gid, g2.gid as gid_ground, g1.x as x, g1.y as y, g2.z as z, g1.z - g2.z as height, g1.the_geom as geometry
 FROM veg As g1, ground As g2   
 WHERE g1.gid <> g2.gid AND ST_DWithin(g1.the_geom, g2.the_geom, 3.5)   
 ORDER BY g1.gid, ST_Distance(g1.the_geom,g2.the_geom);
This is the total query time comparing 167,000 to 222,000 points using ST_DWithin to do a nearest neighbor search with different search distances.

				

Further optimization of the PostGIS LiDAR Vegetation Height Query

There’s much to be said for knowing your data in order to best optimize the analysis of it.  Beyond all other bits of cleverness, having a functional understanding of your problem is the first step toward conceiving an intelligent and efficient solution.

One thing that I didn’t do two posts ago was to spend any time deciding how far out the search for nearby points with ST_DWithin would be.  So I played around a bit visualizing the data to address this question.  Here in green are the ground points shown in Quantum GIS.

And below they are displayed as points two feet in diameter (kind of like buffering them, but here, it’s just taking advantage of qGIS’ display properties).

They almost converge.  Going up to 3 feet, we can conclude that the spacing of this LiDAR dataset is about 3 feet.  These points are over a flat residential area (the rectangular gaps are buildings).

and below we have the points over a cliff/slump area.  We can see some gaps in this extreme topography.  So I’ll conclude that if I use a search area of 3.5 feet, I should be able to more efficiently perform my nearest neighbor search, and find the heights of the vegetation and buildings relative to the closest ground point.

SELECT DISTINCT ON(g1.gid)  g1.gid as gid, g2.gid as gid_ground, g1.x as x, g1.y as y, g2.z as z, g1.z - g2.z as height, g1.the_geom as geometry
 FROM veg As g1, ground As g2   
 WHERE g1.gid <> g2.gid AND ST_DWithin(g1.the_geom, g2.the_geom, 3.5)   
 ORDER BY g1.gid, ST_Distance(g1.the_geom,g2.the_geom);
So, here is the vegetation height shown from lightest green (shortest) to darkest green (tallest):

We can throw in buildings similarly colored in gray (in this case we extend the search window to 30 feet to ensure we find a ground point nearby):

And a closer look over my house.

I know I wrote down some numbers for how much faster this is than searching 15 feet, but I can't remember where I put them.  Needless to say, this is the most important optimization.