Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Posts Tagged ‘SQL’

Plugin-free QGIS TMS tiles via GDAL

Posted by smathermather on July 31, 2014

Want to load your favorite tiles into QGIS? How about a plugin-free QGIS TMS tiles via GDAL:

http://www.3liz.com/blog/rldhont/index.php?post/2012/07/17/OpenStreetMap-Tiles-in-QGIS

Really awesome… .

Needs but one change: epsg:900913 should be epsg:3857 or QGIS (GDAL?) throws an error. Presumably you could also define epsg:900913 in some config file, but barring that create an XML file as follows, and load as a raster in QGIS:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://maps1.clemetparks.com/tilestache/tilestache.cgi/basemap/${z}/${x}/${y}.jpg</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>

Now I can use my pretty tilestache maps _everywhere!_:
Image of coyote points overlayed on tilestache rendered hillshade map

Screenshot of map from postgis

edit: forgot the hattip for the lead on this, Thomas gratier: @ThomasG77

Posted in Analysis, Database, PostGIS, PostgreSQL, QGIS, SQL, TileStache | Tagged: , , , , | 1 Comment »

Cleaning animal tracking data — throwing away extra points

Posted by smathermather on July 31, 2014

Much the problem of the modern era– too much data, uneven data, and yet, should we keep it all?

Here’s the problem space: attach GPS collar to a coyote, send that data home, and you have a recipe for gleaning a lot of information about the movement of that animal across the landscape. In order to maximize the data collected while also maximizing the battery life of said device, sometimes (according to a preset scheme) you might sample every hour or every 15 minutes, and then switch back to once a day.

When you get the data back in, you get something like this:

Map of raw coyote data-- all points, showing all available collar GPS points

Already, we see some pretty interesting patterns. The first thing that I notice is clusters of activity. Are these clusters related to our uneven sampling, sometimes every 15 minutes, sometimes once or twice a day? Or are those really important clusters in the overall activity of the animal?

The next thing I notice is how perfectly the range of movement of the coyote is bounded by expressways to the west and north, and by the river to the east. Those seem to be pretty impermiable boundaries.

So, as does a man who only has a hammer, we clean the data with PostGIS. In this case, it’s an ideal task. First, metacode:

Task one: union / dissolve our points by day — for some days this will give us a single point, for other days a multi-point for the day.

Task two: find the centroid of our union. This will grab a centroid for each clump of points (a weighted spatial mean, if you will) and return that single point per day. There are some assumptions here that might not bear out under scrutiny, but it’s not a bad first approximation.

Now, code:

CREATE TABLE centroids AS
-- union / dissolve our points by day
WITH clusterp AS (
	SELECT ST_Union(geom) AS geom, gmt_date FROM nc_clipped_data_3734
		GROUP BY gmt_date
),
-- find the centroid of our union
centroids AS (
	SELECT ST_Centroid(geom) AS geom, gmt_date FROM clusterp
)
SELECT * FROM centroids;

Boom! One record per day:

Map showing much smaller home range, one record per day clustered around a single core zone

A different pattern now emerges. We can’t see the hard boundaries of the use area, but now the core portion of the home range can be seen.

Upcoming (guest) blog post on the use of these data in home range calculations in R. Stay tuned.

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

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

Editing in progress

Posted by smathermather on January 9, 2014

Edit:

Free book to the first to tell me what the is name of the children’s book on the armchair in the foreground:

 

Book editing in progress… .

image

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

More cutting room floor stuff… — 3D pyramid maker

Posted by smathermather on December 14, 2013

A little more leftover code from the PostGIS Cookbook— a little tool for making pyramids for SFCGAL enabled PostGIS databases.

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

PostGIS Pyramid

Edit:  Actually, I think this would work without SFCGAL, but what would be the point… .

Edit 2:  Let’s embed some code:

 

-- Function creates an inverted pyramid
-- Function takes as input an origin point, the size base in x and y, and the height of the pyramid
-- Usage example: SELECT 1 pyramid_maker(geom, 2, 2, 1)  AS the_geom;
CREATE OR REPLACE FUNCTION pyramidMaker(origin geometry, basex numeric, basey numeric, height numeric)
  RETURNS geometry AS
$BODY$

-- Create points as the base of the pyramid-- which perversly is in the air... .
WITH basePoints AS
	(
	SELECT ST_Translate(origin, -0.5 * basex, 0.5 * basey) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, 0.5 * basex, 0.5 * basey) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, 0.5 * basex, -0.5 * basey) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, -0.5 * basex, -0.5 * basey) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, -0.5 * basex, 0.5 * basey) AS the_geom
	),
-- Make the points into a line so we can convert to polygon
basePointsC AS
	(
	SELECT ST_MakeLine(the_geom) AS the_geom FROM basePoints
	),	
baseBox AS
	(
	SELECT ST_MakePolygon(ST_Force3DZ(the_geom)) AS the_geom FROM basePointsC
	),
-- move base of pyramid vertically the input height
base AS
	(
	SELECT ST_Translate(the_geom, 0, 0, height) AS the_geom FROM baseBox
	),
-- Now we construct the triangles, ensuring that we digitize them counterclockwise for validity
triOnePoints AS
	(
	SELECT origin AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, 0.5 * basex, 0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, -0.5 * basex, 0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT origin AS the_geom
	),
triOneAngle AS
	(
	SELECT ST_MakePolygon(ST_MakeLine(the_geom)) the_geom FROM triOnePoints
	),
triTwoPoints AS
	(
	SELECT origin AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, 0.5 * basex, -0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, 0.5 * basex, 0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT origin AS the_geom
	),
triTwoAngle AS
	(
	SELECT ST_MakePolygon(ST_MakeLine(the_geom)) the_geom FROM triTwoPoints
	),
triThreePoints AS
	(
	SELECT origin AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, -0.5 * basex, -0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, 0.5 * basex, -0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT origin AS the_geom
	),
triThreeAngle AS
	(
	SELECT ST_MakePolygon(ST_MakeLine(the_geom)) the_geom FROM triThreePoints
	),
triFourPoints AS
	(
	SELECT origin AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, -0.5 * basex, 0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT ST_Translate(origin, -0.5 * basex, -0.5 * basey, height) AS the_geom
		UNION ALL
	SELECT origin AS the_geom
	),
triFourAngle AS
	(
	SELECT ST_MakePolygon(ST_MakeLine(the_geom)) the_geom FROM triFourPoints
	),
-- Assemble the whole thing into a pyramid
pyramid AS
	(
	SELECT the_geom FROM triOneAngle
		UNION ALL
	SELECT the_geom FROM triTwoAngle
		UNION ALL
	SELECT the_geom FROM triThreeAngle
		UNION ALL
	SELECT the_geom FROM triFourAngle
		UNION ALL
	SELECT the_geom FROM base
	),
-- format the pyramid temporarily as a 3D multipolygon
pyramidMulti AS
	(
	SELECT ST_Multi(St_Collect(the_geom)) AS the_geom FROM pyramid
	),
-- convert to text in order to manipulate with a find/replace into a polyhedralsurface
-- and then convert back to binary
textPyramid AS
	(
	SELECT ST_AsText(the_geom) AS textpyramid FROM pyramidMulti
	),
textBuildSurface AS
	(
	SELECT ST_GeomFromText(replace(textpyramid, 'MULTIPOLYGON', 'POLYHEDRALSURFACE')) AS the_geom FROM textPyramid
	)

SELECT the_geom FROM textBuildSurface

;

$BODY$
  LANGUAGE sql VOLATILE
  COST 100;
ALTER FUNCTION pyramid_maker(geometry, numeric, numeric, numeric)
  OWNER TO me;

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

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;

Posted in 3D, Database, PostGIS, PostgreSQL, SFCGAL, SQL | Tagged: , , , , | 4 Comments »

Multi-line Comments in SQL for PostgreSQL

Posted by smathermather on November 3, 2012

Who knew? Not I! I shouldn’t confess, but just discovered multi-line comments in PostgreSQL:

-- This is a single line comment in SQL (I knew this already...)
/* Multi-line quotes use C
syntax.  Simple as that.  Surprised I hadn't googled (well duck duck goed) this before! */

Sadly, WordPress syntax highlighting doesn’t recognize this… . I think I can muddle through without it though.

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

Proper (ab)use of a database, contour interpolation using #postgresql #postgis #arcgis

Posted by smathermather on August 6, 2012

Anyone who has been following along at home knows I don’t think much like a DBA.  Sometimes that’s good; mostly it’s probably bad.  In this post, I hope it will be interesting.

The problem of the day is how to take engineering contours derived from breaklines, a lidar point cloud, and all the lot, and do a good job interpolating that to a DEM.  This is an interesting problem, as it’s another one of those problems with so many sub-optimal solutions.  In part, it’s an impossible problem– going from a surface model to contours is inherently a loss of information, so going the other way is difficult to do believably, both with respect to believable smoothness and the requisite feature retention.  This is a little like the pugilistic pie puppets problem: ‘flaky vs. tender’ of making good pie crust.  Doing one well seems to preclude doing the other well.

Now, Chris Hormann, of PovRay fame has a solution, there are solutions in ArcGIS, and in GRASS, and plenty of other attempts.  Here’s a decent list, which I haven’t checked for timeliness:  http://vterrain.org/Elevation/contour.html.  Since we have a ArcGIS license, it handles the dataset I have OK, the results checked out OK, this is what we chose to use for now.  I tried GRASS, but it was too slow for the amount of data being processed.  The command in ArcGIS ends up looking something like this:


TopoToRaster_sa "N2260635.shp elevation Contour;N2260640.shp elevation Contour;N2260645.shp elevation Contour;N2265633.shp elevation Contour;N2265640.shp elevation Contour;N2265645.shp elevation Contour;N2270635.shp elevation Contour;N2270640.shp elevation Contour;N2270645.shp elevation Contour;N2275650.shp elevation Contour" "C:\Temp\TopoToR_N2269.tif" 2 "2265031 645070 2270065 639945" 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # # #

And thus we go from this:

To this:

Now to do it 623 more times… .  And in walks PostgreSQL for some abuse as a spatial engine with string manipulation to boot.

Each tile of contours I run through this approach needs the adjacent 8 tiles of contours loaded in order to avoid dreaded edge effects:

Except that sometimes, there aren’t 8 tiles available, like on the edge of large water bodies:

So, we can’t just use a heuristic that loads the 8 adjacent tiles– we need to enumerate these tiles.  An ST_Touches should do the trick here:

SELECT p.tile, a.tile AS tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom);

But, that does not arrange our table very nicely:

I want something that looks like this:

N2225655, N2225660, N2220650 …

N2225660, N2225665, N220665 … etc.

Thank goodness for Postgres’ String_agg, and moreover for Leo Hsu and Regina Obe’s blog.

Adapting their entry, we can reorganize those entries like this:

SELECT
STRING_AGG(a.tile, ',' ORDER BY a.tile)
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile;

Ah, now that’s better:

And if we want, we can use an arbitrary delimiter, and build our entire command in one swell foop:

SELECT p.tile,
'TopoToRaster_sa "'
|| 'W:\contours\cuy_contours\shapefiles\' || p.tile || '.shp elevation Contour;'
|| 'W:\contours\cuy_contours\shapefiles\'
|| STRING_AGG(a.tile, '.shp elevation Contour;W:\contours\cuy_contours\shapefiles\' ORDER BY a.tile)
|| '.shp elevation Contour;" '
|| '"v:\svm\dem\cuyahoga\' || p.tile || '.tif" 2'
|| ' "' || replace(replace(replace(ST_Extent(p.geom)::text, 'BOX(', ''), ')', ''), ',', ' ') || '"'
|| ' 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # '
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile

Ok.  I have 623 more tiles to do.  Until next time… .

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

Going deeper into web cartography: future=past? (and Swiss cartographic genius)

Posted by smathermather on May 12, 2012

My favorite cartography book is Eduard Imhof’s Cartographic Relief Presentation.  A few years back I picked this book up (translated to English) from ESRI press for $75 if memory serves me.  Now it can be gotten for much cheaper.

Not created in GeoServer (nor TileMill– from Imhof’s book)

Imhof spends a lot of time on feature simplification and separation, a problem which keeps me up at night.  For example, if you have a lot of overlapping and/or touching contours, what are the local distortions and simplifications that can be done to enhance map readability?

Cartographic representation choices with overlapping and tightly packed contours, also from Cartographic Relief Presentation

This problem applies to other overlapping features, such as a road following a river, where at a given scale one might not distinguish between the two, so we exaggerate their differences, in this case, while maintaining their basic topology.

I’ve played with this problem before, but now I’m considering a new approach.  I pinged Martin Davis of JTS fame for some advice.  He suggested using force directed layout to solve this problem, much as he does in JTS Test Builder for magnifying/investigating topology.

I don’t know yet where or how (if) I’ll implement this, but it’s an interesting an potentially solvable problem.  Now where to work on this in the stack– in PostGIS, because I know it best, in the renderer (GeoServer) where I’ll be swimming in Java of a level I don’t have a chance, or in the browser/client, where there are already some examples written in Javascript… .

Posted in Analysis, Cartography, Database, PostGIS, PostgreSQL, SQL, Trail Curation, Trails | Tagged: , , , , , , , | Leave a Comment »

PostGIS for Dessert: Sketching shapes in #PostGIS– compass roses revisited

Posted by smathermather on May 10, 2012

For one of our applications, we need 8-point compass roses placed at each of our points, as well as a circle 40 meters in diameter as well as one 140 meters in diameter.  We did a bit of work with this a while back in GeoServer using SLDs.  Now we’d like to refine it, and implementing this in an SLD is beyond my skills.

So, we move to the back end.  The following is a little function to construct the roses for us, which we’ll then display though GeoServer.  Just an FYI, this is running on a PostGIS 1.3 database, and I didn’t (at the time anyway) grok the use of ST_RotateZ, so there are some explicit position calls, rather than just creating the crosshairs and rotating it:


CREATE OR REPLACE FUNCTION pie (
  geom Geometry, filling numeric, crust numeric
 )
 RETURNS Geometry
 AS $$
DECLARE

  slices Geometry;

BEGIN

	--slices = ST_MakeLine(geom, ST_Translate(geom, dough, pan));
	slices =

	ST_Union(
		ST_Collect(
			ST_Collect(
				ST_Collect(
					ST_MakeLine(ST_Translate(geom, filling, 0), ST_Translate(geom, crust, 0))
					,ST_MakeLine(ST_Translate(geom, 0, filling), ST_Translate(geom, 0, pi()/4 * crust))
				)
				,ST_Collect(
					ST_MakeLine(ST_Translate(geom, -filling, 0), ST_Translate(geom, pi()/4 * -crust, 0))
					,ST_MakeLine(ST_Translate(geom, 0, -filling), ST_Translate(geom, 0, pi()/4 * -crust))
				)
			)
			,ST_Collect(
				ST_Collect(
					ST_MakeLine(ST_Translate(geom, filling * pi()/4, filling * pi()/4), ST_Translate(geom, crust * pi()/4, crust * pi()/4))
					,ST_MakeLine(ST_Translate(geom, -filling * pi()/4, filling * pi()/4), ST_Translate(geom, -crust * pi()/4, crust * pi()/4))
				)
				,ST_Collect(
					ST_MakeLine(ST_Translate(geom, filling * pi()/4, -filling * pi()/4), ST_Translate(geom, crust * pi()/4, -crust * pi()/4))
					,ST_MakeLine(ST_Translate(geom, -filling * pi()/4, -filling * pi()/4), ST_Translate(geom, -crust * pi()/4, -crust * pi()/4))
				)
			)
		)
		,ST_Collect(
			ST_ExteriorRing(ST_Buffer(geom, filling, 32))
			,ST_ExteriorRing(ST_Buffer(geom, crust, 32))			
		)
	)
	
		;
	
RETURN slices;
END;

$$ LANGUAGE plpgsql;

So now to test it with a single point:


DROP TABLE IF EXISTS test_pie;

CREATE TABLE test_pie AS

                SELECT pie(ST_MakePoint(0,0), 40, 140);

And view in uDig, or application of choice:

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

Building simple clients for MapFish — cURL as a client

Posted by smathermather on April 25, 2012

I have two previous posts on using MapFish (in this case, the GeoServer version) to allow for printing to hi-resolution PDF maps from the browser.  Here we use a command-line browser (cURL) to post our json to the MapFish service in order to retrieve our PDF.

I did not keep any notes from before on making json posts to the MapFish server as a means by which to test any manual configuration of the json file, so I had to rediscover this approach from pages like this.

The “@” sign below is so that curl knows I’m feeding it a file instead of the actual json to post:


curl -i -H "Accept: application/json" -H "Content-Type: application/json" -X POST --data @mapfish_landscape.json http://localhost:8080/geoserver/pdf/create.json

Posted in Database, GeoExt, GeoExt, GeoServer, MapFish, PostGIS, PostgreSQL, SQL | Tagged: , , , , , , , , , | 3 Comments »