Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for December, 2013

2013 in review

Posted by smathermather on December 30, 2013

The WordPress.com stats helper monkeys prepared a 2013 annual report for this blog.

Here’s an excerpt:

The concert hall at the Sydney Opera House holds 2,700 people. This blog was viewed about 38,000 times in 2013. If it were a concert at Sydney Opera House, it would take about 14 sold-out performances for that many people to see it.

Click here to see the complete report.

Posted in Fluff | Tagged: | Leave a Comment »

Drone Pointilism

Posted by smathermather on December 23, 2013

Just an image today.

image of UAS derived points using structure from motion.

Image of UAS derived points using structure from motion.

Posted in 3D, Bundler, Drone, Image Processing, Optics, Photogrammetry, PMVS, PostGIS, QGIS, UAS | Tagged: , , , | Leave a Comment »

Voronoi in PostGIS

Posted by smathermather on December 21, 2013

PostGIS has for a time had ST_DelaunayTriangles for Delaunay Triangulation, and since 2.1 apparently can even natively create a 2.5D TIN using said function, which is pretty cool. I think with SFCGAL, we will eventually have true 3D TINs as well.

Overlay of Delaunay and Voronoi

We’re still waiting on native Vororoi support for PostGIS though. According to http://trac.osgeo.org/postgis/wiki/UsersWikiPostgreSQLPostGIS

“GEOS 3.5 is still in development but will have Voronoi which will be a feature in 2.2 only if you have GEOS 3.5.”

Vishal Tiwari through a Google of Summer of Code under the mentorship of Sandro Santilli completed the port of Voronoi to GEOS from JTS Topology Suite. Now we just need to wait for PostGIS 2.2….

In the mean time,

http://geogeek.garnix.org/2012/04/faster-voronoi-diagrams-in-postgis.html

One caveat– this python function doesn’t always provide just Voronoi but some artifact polygons.

Voronoi polygons with some artifacts

Once we have a table with Voronoi, we can filter out just the true Voronoi cells by counting the number of original points we find within them, and only return the polygons which contain a single point:

CREATE TABLE true_voronoi AS
WITH tvoronoi AS (
	SELECT COUNT(*), v.geom
		FROM voronoi v, input_points p
		WHERE ST_Intersects(v.geom, p.geom)
		GROUP BY v.geom
		)
SELECT the_geom FROM tvoronoi WHERE count = 1;

voronoi without extra polygons

But it’s still not a perfect solution. I can’t wait for PostGIS 2.2….

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

What’s in a name? That which we call a GIS By any other name would smell as sweet;

Posted by smathermather on December 20, 2013

What’s in a name? That which we call a GIS
By any other name would smell as sweet;

But would it. I’m at a place where I could name my working group. Thoughts? (we could just call ourselves the GeoHipsters…).

GIS, Geospatial, Mapping…

Posted in Other | Tagged: | 2 Comments »

They’ve gone to plaid. (LeafletJS and CartoDB)

Posted by smathermather on December 18, 2013

Enjoyed setting the errorTileUrl flag in a leafletjs map today.

errorTileUrl: 'https://encrypted-tbn2.gstatic.com/images?q=tbn:ANd9GcStyzXdR5V4U_BF_DOOBAl3eEJGeNw54O40_Gxj5nXRMc-49Ega0A'

Image of map with plaid background

A small homage:

Posted in Leaflet, Other | Tagged: , , | 2 Comments »

2.5D TINs in PostGIS

Posted by smathermather on December 18, 2013

(edited: changed TIN to TIN Z)
(edited again — function already exists as a flag in ST_DelaunayTriangles… .)

Uh, I wrote a function for nothin’…
As Regina points out in the commments, this functionality was rolled into 2.1 with a flag. Her example:

SELECT
    ST_AsText(ST_DelaunayTriangles(
       ST_Union(ST_GeomFromText(
       'POLYGON((175 150, 20 40, 50 60, 125 100, 175 150))'), 
       ST_Buffer(ST_GeomFromText(‘POINT(110 170)’), 20) ),0.001,2) )
    AS dtriag;

For the record, my code is about 2% slower on a dataset that takes ~50minutes to triangulate .

——————————————
Original Post
——————————————

Ever since ST_Delaunay made its way into PostGIS with GEOS 3.4, there’s been the promise of TINs from geometries in PostGIS. The promise is for 2.5D TINs.  What are these, and why aren’t they 3D? It’s all about overhangs– imagine the edge of a building.  A 3D TIN can handle that building edge with a true vertical wall, or even overhang.  A 2.5D TIN is like a 2.5D raster– no overlaps allowed.

With those limitations in mind, you can have TINs today if you want:
https://github.com/smathermather/postgis-etc/edit/master/3D/AsTin.sql

-- A small function to convert ST_Delaunay (requires GEOS 3.4) to a 2.5D Tin
-- Uses the hackerish approach of converting to text, and doing string replacement
--- for format conversion.

-- A Delaunay triangulation does not by itself formally make a TIN.  To do this,
-- we will require some modification to the output format.  If we perform an ST_AsText
-- on our new geometries we will see that they are POLYGON Z's wrapped in a
-- GEOMETRYCOLLECTION.

-- Thus there are two ways in which this is not a TIN-- the internal POLYGON Z are
-- implicitly triangles, but explicitly POLYGON Zs.  In addition, the external wrapper
-- for the collection of triangles is a GEOMETRYCOLLECTION, not a TIN.
-- Once that we have this geometry in text form, the easiest way to fix this is
-- with string replacement to fix these two things, then convert back to binary

-- We'll go from e.g. this:
-- GEOMETRYCOLLECTION Z (POLYGON Z ((-1.14864 0.61002 -2.00108,-0.433998 -0.0288903 -2.13766, ...
-- to this:
-- TIN ( ((-1.14864 0.61002 -2.00108,-0.433998 -0.0288903 -2.13766, ...

CREATE OR REPLACE FUNCTION AsTIN(geometry)
  RETURNS geometry AS
$BODY$

WITH dt AS
(
SELECT ST_AsText(ST_DelaunayTriangles(ST_Collect($1))) AS atext
),
replacedt AS
(
-- Remove polygon z designation, as TINs don't require it.
SELECT replace(atext, 'POLYGON Z', '') as ttext
  FROM dt
),
replacegc AS
(
-- change leading declaration to TIN
SELECT replace(ttext, 'GEOMETRYCOLLECTION Z', 'TIN Z') AS tintext
  from replacedt
),
tingeom AS
(
-- Aaaand convert back to binary.  Voila!
SELECT ST_GeomFromEWKT(tintext) AS the_geom FROM replacegc
)

SELECT the_geom FROM tingeom

$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

Here’s an example TIN derived from UAS (non-weaponized drone) imagery:
Image showing UAS derived TIN

Posted in 3D, Analysis, Database, Drone, Other, Photogrammetry, PostGIS, PostgreSQL, SQL, UAS | Tagged: , , , , , | 6 Comments »

UAS (drone) Footprint Geometries Calculated in PostGIS with SFCGAL — for real this time

Posted by smathermather on December 15, 2013

In my earlier post, I made a claim that SFCGAL was used in this figure:

Figure showing overlapping viewing footprints from images from UAS flight.

It dawned on my afterwards, while I was using 3D, I hadn’t actually employed any of the analysis goodies that come with SFCGAL.  Well, here it is– a footprint as calculated using the view angles and a real terrain model:

Map showing TIN of intersection of digital terrain model and viewing cone from UAS.

Here it is compared with the initial calculation:

Comparison of previous figure with original less correct estimate.

Posted in 3D, Analysis, Database, Drone, Image Processing, Optics, Other, Photogrammetry, PostGIS, PostgreSQL, QGIS, SFCGAL, SQL, UAS | Tagged: , , , , , | 3 Comments »

UAS (drone) Footprint Geometries Calculated in PostGIS — Individual Footprint

Posted by smathermather on December 15, 2013

UAS (drone) Footprint Geometries Calculated in PostGIS (viewed in QGIS nightly), taking into account relative elevation, bearing, pitch, and roll, this time just one:

Nadir view of viewing pyramid for individual drone  image

Posted in 3D, Analysis, Database, Drone, Image Processing, Optics, Other, PostGIS, PostgreSQL, QGIS, SFCGAL, SQL, UAS | Tagged: , , , | Leave a Comment »

UAS (drone) Footprint Geometries Calculated in PostGIS with SFCGAL

Posted by smathermather on December 15, 2013

UAS (drone) Footprint Geometries Calculated in PostGIS (viewed in QGIS nightly), taking into account relative elevation, bearing, pitch, and roll:

Figure showing overlapping viewing footprints from images from UAS flight.

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

Aaargh! No: geometry ST_RotateX(geometry geomA, float rotRadians, geometry pointOrigin)

Posted by smathermather on December 14, 2013

Edit– Code may be flawed, testing—– testing. Please wait to use… .

In PostGIS, ST_RotateZ has a couple forms: a version that rotates around 0,0 and a version that rotates around a point of origin:

geometry ST_Rotate(geometry geomA, float rotRadians, geometry pointOrigin)

ST_RotateX and ST_RotateY have had no equivalents– until now. These equivalents are dumb versions– they use a transformation/rotation/reverse-transformation to do their magic, which is maybe not as efficient as using ST_Affine, but I’m not smart enough for ST_Affine. Not today anyway. Repo here: https://github.com/smathermather/postgis-etc/tree/master/3D

geometry ST_RotateX(geometry geomA, float rotRadians, geometry pointOrigin)

ST_RotateX version:

-- Function: st_rotatex(geometry, double precision, geometry)
CREATE OR REPLACE FUNCTION ST_RotateX(geomA geometry, rotRadians double precision, pointOrigin geometry)
  RETURNS geometry AS
$BODY$

----- Transform geometry to nullsville (0,0,0) so rotRadians will take place around the pointOrigin
WITH transformed AS (
    SELECT ST_Translate(geomA, -1 * ST_X(pointOrigin), -1 * ST_Y(pointOrigin), -1 * ST_Z(pointOrigin)) AS the_geom
    ),
----- Rotate in place
rotated AS (
    SELECT ST_RotateX(the_geom, rotRadians) AS the_geom FROM transformed
    ),
----- Translate back home
rotTrans AS (
    SELECT ST_Translate(the_geom, ST_X(pointOrigin), ST_Y(pointOrigin), ST_Z(pointOrigin)) AS the_geom
	FROM rotated
    )
----- profit
SELECT the_geom from rotTrans

;

$BODY$
  LANGUAGE sql VOLATILE
  COST 100;
geometry ST_RotateY(geometry geomA, float rotRadians, geometry pointOrigin)

ST_RotateY version:

-- Function: ST_RotateY(geometry, double precision, geometry)
CREATE OR REPLACE FUNCTION ST_RotateY(geomA geometry, rotRadians double precision, pointOrigin geometry)
  RETURNS geometry AS
$BODY$

----- Transform geometry to nullsville (0,0,0) so rotRadians will take place around the pointOrigin
WITH transformed AS (
    SELECT ST_Translate(geomA, -1 * ST_X(pointOrigin), -1 * ST_Y(pointOrigin), -1 * ST_Z(pointOrigin)) AS the_geom
    ),
----- Rotate in place
rotated AS (
    SELECT ST_RotateY(the_geom, rotRadians) AS the_geom FROM transformed
    ),
----- Translate back home
rotTrans AS (
    SELECT ST_Translate(the_geom, ST_X(pointOrigin), ST_Y(pointOrigin), ST_Z(pointOrigin)) AS the_geom
	FROM rotated
    )
----- profit
SELECT the_geom from rotTrans

;

$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

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