Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

More cutting room floor stuff… .

Posted by smathermather on December 14, 2013

ST_3DIntersection performs volumetric intersections for us in PostGIS, if we have SFCGAL enabled for our PostGIS back end.  Much like a normal PostGIS intersection, this is the mathematical definition of an interesection, so it returns the volumetric portion of the intersection, plus 3D linestrings and 3D points and other bits and pieces that qualify for the intersection.  As a little patch, I wrote a quick and dirty function for just extracting the volumetric portion.  It’s available here:

https://github.com/smathermather/postgis-etc/blob/master/3D/volumetricIntersection.sql

It’s a little cludgy, as not all the bits and pieces we’re used to on the 2D side are built yet, but it works!

One of these days I’ll get horao working on my machine… .

Edit: I’ll have more on volumetric intersection to come, but in the meantime, this tutorial on PovRay:
http://www.cs.tut.fi/~tgraf/harjoitustyot/tutorial/tutorial1.6.html does a good job explaining.

For example, if we intersect two overlapping spheres:
Figure of overlapping spheres
then we get an output thus (a volumetric intersection of two spheres):
Figure of the intersection result of two spheres

edit2, relevant code:

CREATE OR REPLACE FUNCTION volumetricIntersection(geom1 geometry, geom2 geometry)
-- volumetric intersection takes an the input of two 3D geometries
  RETURNS geometry AS
$BODY$

WITH 
intersected AS (
-- first we perform an intersection. This in most cases will return a TIN plus 3D linstrings
-- and other messy pieces we don't need.
	SELECT ST_3DIntersection(geom1, geom2) AS the_geom
	),
-- we use ST_Dump to dump these out to their requisite parts
-- (no, ST_CollectionExtract will not work here-- it only handles
-- points, lines, and polygons, not triangles and tins
dumped AS (
	SELECT (ST_Dump(the_geom)).geom AS the_geom FROM intersected
	),
-- Now we filter for triangle and collect them together
triangles AS (
	SELECT ST_Collect(the_geom) AS the_geom FROM dumped
		WHERE ST_GeometryType(the_geom) ='ST_Triangle'
	),
-- next as a venerable hack, we'll convert to text
triangleText AS (
	SELECT ST_AsText(the_geom) AS triText FROM triangles
	),
-- and replace words in the text in order to "convert" from a collection
-- of triangles to a TIN
replaceTriangle AS (
	SELECT replace(triText, 'TRIANGLE Z ', '') AS rTri FROM triangleText
	),
tinText AS (
	SELECT replace(rTri, 'GEOMETRYCOLLECTION', 'TIN') AS tt FROM replaceTriangle
	),
-- now we convert back to a binary tin, and give it back to the user
tin AS (
	SELECT ST_GeomFromText(tt) AS the_geom FROM tinText
	)
	
SELECT the_geom FROM tin;

$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

4 Responses to “More cutting room floor stuff… .”

  1. mapbaker said

    this one… I just don’t get…

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: