Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for the ‘SQL’ Category

Using PostGIS for Hydrologic Modeling (reblog)

Posted by smathermather on June 17, 2016

The Problem We have to filter out the roads and ditches without removing streams that cross roads or follow them closely. I’m going to use PostGIS to find the intersection of the streams lines data with a buffered roads polygon. If the intersected line is less than 50% of the length of the stream line, […]

via Filtering Roads from Extracted Streams Data — GeoKota

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

Using foreign data wrapper to use PostGIS with SQLServer

Posted by smathermather on May 29, 2016

Here was the problem that needed solved last week (we have a few similar problems in upcoming projects, so this was an exciting thing to try): we needed to use PostGIS to access data in a SQLServer database. The SQLServer database backs the web site in question, the underlying content management system, etc., so no– removing SQLServer isn’t really an option at this stage. Obviously PostGIS is a requirement too… .

Before I go further, I used tds_fdw as the foreign data wrapper. The limitations here are as follows: it is a read-only connection, and only works with Sybase and SQLServer, as it uses tabular data stream protocol for communicating between client and server. This is not as generic a solution as we can use. Next time I’ll try ogr_fdw which is more generic as it can connect with other databases and other data types. Another advantage to ogr_fdw is we can use IMPORT FOREIGN SCHEMA. Regina Obe and Leo Hsu warn though to limit this to 150 tables or so for performance reasons.

With the limitations listed above, this is how we built the thang:

DROP SERVER beach_fdw CASCADE;

-- Create the server connection with FOREIGN DATA WRAPPER
CREATE SERVER beach_fdw
FOREIGN DATA WRAPPER tds_fdw
OPTIONS (servername 'name_or_ip', port '1433', database 'non-spatial-database', tds_version '7.1', msg_handler 'notice');

-- We map the postgres user to the user that can read the table we're interested in
CREATE USER MAPPING FOR postgres
SERVER beach_fdw
OPTIONS (username 'user', password 'password');

-- Create the actual foreign table connection
CREATE FOREIGN TABLE beach_closures (
AutoNumber int NOT NULL,
Title varchar NOT NULL,
StartDate timestamp without time zone NOT NULL,
WaterQuality varchar NOT NULL,
Latitude varchar NOT NULL,
Longitude varchar NOT NULL,
BeachStatus varchar NOT NULL,
ClosureColor varchar NOT NULL)
SERVER beach_fdw
OPTIONS (schema_name 'schema_name', table_name 'vw_CMPBeachClosures');

-- Now we create a spatial view using our longitude and latitude
CREATE VIEW v_beach_closures AS
SELECT
AutoNumber, Title, StartDate, WaterQuality, Latitude,
Longitude, BeachStatus, ClosureColor, ST_SetSRID(ST_MakePoint(Longitude::numeric, Latitude::numeric), 4326)	AS geom
FROM beach_closures;

Voila! A nice little PostGIS enabled view of a SQLServer view or table!

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

Point Clouds – the (re)start of a journey

Posted by smathermather on January 19, 2016

If you have followed this blog for a while, you will notice a continual returning to and refinement of ideas and topics. That’s how the blog started, and this role it has served, as a touch stone in my exploration of topics is critical to how I use it. I hope it is useful to you too, as a reader.

So, let’s talk about point clouds again. Point clouds are like ordinary geographic points — they have X, Y, and Z (typically), but there are a whole lot more of them than ordinary. How much more? Well, instead of thinking in the 100s of thousands, we think more along the lines of 100 of millions or billions of points.

The other thing we need to recognize about point clouds is that they may be n-dimensional — a LiDAR point cloud may have X, Y, Z, intensity, return number, class, etc..

Good Free and Open Source Software (FOSS) tools for dealing with point clouds include a spectrum of projects, including PDAL (Point Data Abstraction Library), and the Pointcloud extension for PostgreSQL (with PostGIS integration, of course).

Previously, I was attempting to process large amounts of LiDAR point clouds in order to model bird habitat. Now that I have the requisite storage, and PDAL has a height calculator, I am ready to dive back into this.

In the meantime, there are some cool things to check in this space. Let’s start with this presentation by Michael Smith on PDAL from FOSS4G in Seoul last year.

Screen Shot 2016-01-19 at 8.32.03 PM

Want to see some of the incredible things you can do in the browser, go no further than here:

Screen Shot 2016-01-19 at 8.34.52 PM.png

This goes deeper than simple vector tile approaches and employs level of detail optimizations that are necessary for this large of a use case. More on this later

postscript 1/20/2016:

Looks like I outed the above project… . Apologies to https://hobu.co for this. For the record and from their site:

Hobu is a team of three software engineers located in Iowa City, Iowa and Cambridge, Iowa. We have been on the forefront of open source LiDAR software for over five years, and we have been building open source GIS software for even longer. We provide open source lidar software systems development, open source GIS software development in addition to contract research and development, systems design, evaluation and implementation.

Please contact us about any of the open source software for which we provide support.

We are the primary support organization for the following libraries:

  • PDAL — Point Data Abstraction Library
  • plas.io — WebGL point cloud rendering
  • Greyhound — Point cloud data streaming
  • libLAS — ASPRS LAS read/write

Welp, there we go, an unintentional public launch. Awesome work though… .

Posted in 3D, Database, FLANN, pointcloud, PostGIS, PostgreSQL, SQL | Tagged: , , | 1 Comment »

Generative art in PostGIS

Posted by smathermather on January 14, 2016

Screen Shot 2016-01-14 at 11.46.55 PMScreen Shot 2016-01-14 at 11.45.55 PMScreen Shot 2016-01-14 at 11.45.41 PM

SELECT ST_Buffer(ST_SubDivide(ST_Buffer(
     ST_GeomFromText(
      'POINT(0 0)'
     ), 50000, 100), generate_series), generate_series * 100) AS geom
	FROM generate_series(8,100)

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

Yet another approach to ST_Buffer on geography

Posted by smathermather on January 13, 2016

Another approach to ST_Buffer would be to subdivide the geometries before buffering, and put them all together at the end. ST_SubDivide can do this for us. We can tell it how may vertices we want in each geometry (minimum of 8). Since _ST_BestSRID will try UTM first, we’ll add enough nodes to ensure we always have 8 nodes within the 1,000,000 meter width of a UTM zone by segmentizing at 62,500 meters.

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$

WITH subdivided AS (
 SELECT
 CASE WHEN ST_GeometryType(g1::geometry) = 'ST_Point' THEN g1
 ELSE (SELECT ST_SubDivide(ST_Segmentize(g1::geometry, 62500), 8) AS geom)
 END
),
transformed_local AS (
 SELECT ST_Transform(g1::geometry, _ST_BestSRID(geom)) AS geom FROM subdivided
),
buffered AS (
 SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
 FROM transformed_local
),
transformed_4326 AS (
 SELECT (ST_Dump(ST_Transform(geom, 4326))).geom AS geom FROM buffered
),
checksrid AS (
 SELECT geom FROM transformed_4326
)

SELECT ST_Union(geom)::geography FROM checksrid;
 
$BODY$
LANGUAGE sql VOLATILE;

Aaand, a somewhat unrelated image (just the buffered subdivision of a buffered point, not really an output from our function above. More on all this later.

Screen Shot 2016-01-13 at 6.38.15 PM

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

ST_Buffer on Geography — Iteration

Posted by smathermather on January 12, 2016

In my previous post, Imagico said:

This is nice but still fails in one important aspect – it will create significant inaccuracies once your buffering distance reaches values where buffering in a single local coordinate system is getting problematic which, depending on your precision requirements, will likely happen at distances of a few degrees. A solution to this would require either an iterative split and merge technique of a completely different approach on a very basic level (i.e. attempting truly spherical buffering).

Ok, ok. I kinda forgot about this aspect. The proper way to do this is, as always is to do it at a lower level thank in an SQL function, preferably calculating the buffer on the spheroid, ellipsoid, or geoid. But, ouch. I think we can munge something good enough for now.

We’ll still do the buffer:

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
...

But now that we’ve performed the buffer, we need to determine that the buffered geometry is not too large for our local coordinate system. To do that, we’ll compare the best SRID for our buffered geometry to the original best SRID that we calculated based on the original input geometry.

SELECT
CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography

If the two SRIDs aren’t equal, then we need to use the new best SRID based on buffered geometry and buffer using that new best SRID.

ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, num_seg_quarter_circle), 4326)::geography

Now let’s put it all together:

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$
 
WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT
 CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography 
 ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, num_seg_quarter_circle), 4326)::geography
 END
 FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;


CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
buffer_style_parameters text)
RETURNS geography AS
$BODY$
 
WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, buffer_style_parameters) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT
 CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography 
 ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, buffer_style_parameters), 4326)::geography
 END
 FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;

So, what’s the affect? Pretty significant really.

SELECT ST_Buffer(
	 ST_GeomFromText(
	  'POINT(0 0)'
	 )::geography, 5000000, 50);
Blue line showing correct buffer, dashed red line showing distorted buffer from using wrong local coordinate system.

Blue line showing correct buffer, dashed red line showing distorted buffer from using wrong local coordinate system.

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

Final?: ST_Buffer on Geography, or My love for @NullIsland

Posted by smathermather on January 11, 2016

Thanks to Paul Norman’s reminder in a previous post, we now have all the pieces we need to complete an ST_Buffer function that exposes all the wonderful goodness of buffer style parameters to geography users and also chooses the best local coordinate system automatically.

We use _ST_BestSRID to choose our local coordinate system.

Remember my previous three disadvantages:

  1. It didn’t do this at a low enough level to automatically use the best available local coordinate system. This leaves it to the user to choose the best available local coordinate system
  2. It used a different function name than you otherwise use for buffering.
  3. Finally, it wouldn’t let you use all the other parameters that ST_Buffer exposes through the use of buffer style parameters.

They are now all solved.


CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT geom::geography FROM transformed_4326

$BODY$
LANGUAGE sql VOLATILE;

 

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
buffer_style_parameters text)
RETURNS geography AS
$BODY$

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, buffer_style_parameters) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT geom::geography FROM transformed_4326

$BODY$
LANGUAGE sql VOLATILE;

Let’s try this!

SELECT ST_Buffer(
ST_GeomFromText(
'LINESTRING(5 2.5,0 -2.50,-5 2.5)'
)::geography, 500000, 'join=mitre');
Heart shaped buffer over Null Island.

Heart shaped buffer over Null Island.

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

(poor editing earlier tonight. My apologies.

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

ST_Buffer diving under the covers part 1

Posted by smathermather on January 10, 2016

This will be a quick post. I just looked at how ST_Buffer is written in the PostGIS codebase, and I thought this was interesting. When you use ST_Buffer and specify geometry, buffer distance, and an integer value, this gets converted in the overloaded function into a quad_segs version.


CREATE OR REPLACE FUNCTION ST_Buffer(geometry,float8,integer)

RETURNS geometry

AS $$ SELECT _ST_Buffer($1, $2,

CAST('quad_segs='||CAST($3 AS text) as cstring))

$$

LANGUAGE 'sql' IMMUTABLE STRICT;

P.S. I would have written more, but Sherlock just started.

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

ST_Buffer for those who want really want all buffer features in geography

Posted by smathermather on January 9, 2016

My previous post on improving buffering geography in PostGIS has three disadvantages:

  1. It doesn’t do this at a low enough level to automatically use the best available local coordinate system. This leaves it to the user to choose the best available local coordinate system. I should fix this (but I probably won’t just now…).
  2. It uses a different function name than you otherwise use for buffering. Can we use the same name as ST_Buffer?
  3. Finally, it won’t let you use all the other parameters that ST_Buffer exposes through the use of buffer style parameters.

First, let’s address point 2. With PostgreSQL functions, it is possible to “overload” a function. This means that the same function name, but with a different set of inputs is effectively a different function but with the same name. This is because function + inputs = signature. Since our inputs for our new function are different in number and type from the existing ST_Buffer functions, we can add our function using the same name. The disadvantage is that we will likely lose and have to re-create these functions when we dump and reload the database. (For the advanced user, put it in an alternate schema, and add that to your schema search path so it’s readily available when you want to use it, but you don’t have to go out of your way to make sure that it dumps and reloads.)


CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
 num_seg_quarter_circle integer, local_coord integer)
RETURNS geometry AS
$BODY$
 
WITH transformed AS (
 SELECT ST_Transform(g1::geometry, local_coord) AS geom
),
buffered AS (
 SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
 FROM transformed
),
transformed_4326 AS (
 SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT * FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;

We can also address point 3. Let’s perform a similar overload for buffer style parameters:

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
 buffer_style_parameters text, local_coord integer)
RETURNS geometry AS
$BODY$
 
WITH transformed AS (
 SELECT ST_Transform(g1::geometry, local_coord) AS geom
),
buffered AS (
 SELECT ST_Buffer(geom, radius_of_buffer, buffer_style_parameters) AS geom
 FROM transformed
),
transformed_4326 AS (
 SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT * FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;

Finally, we can use our special buffer style parameters with geography:

SELECT 0 AS id, ST_Buffer(
 ST_GeomFromText(
 'LINESTRING(41 -81,41 -82,42 -82)'
 ), 500000, 'join=mitre mitre_limit=5.0', 3734);

 

Example buffer of line with join=mitre mitre_limit=5.0.

Example buffer of line with join=mitre mitre_limit=5.0.

This is a hack. Daniel’s note in the previous post still stands — it would be better to do this at a lower level within geography, rather than asking the user to know and employ the best possible local coordinate system.

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

ST_Buffer for those who want really round geographic circles

Posted by smathermather on January 8, 2016

A little functionality request from @antoniolocandro on twitter:

Ask and ye might receive. Let’s build a function that will take your input geography, how far you want to buffer (in local coordinates) number of segments you want in a quarter, and what local coordinate system you want to use for your buffer as an EPSG code.

CREATE OR REPLACE FUNCTION zz_buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer, local_coord integer)
RETURNS geometry AS
$BODY$

WITH transformed AS (
SELECT ST_Transform(g1::geometry, local_coord) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT * FROM transformed_4326

$BODY$
LANGUAGE sql VOLATILE
COST 100;

Now let’s use this:

SELECT 0 AS id, zz_buffer(ST_GeomFromText ('POINT(-81 41)', 4326) , 15000, 2, 3734)
Image of 8 segment buffer from geography

Image of 8 segment buffer from geography.

Now, we can test which level of smoothness we like best:

SELECT 0 AS id, zz_buffer(ST_GeomFromText ('POINT(-81 41)', 4326) , 15000, generate_series, 3734)
FROM generate_series(1,50)
Square, Octagon, etc. generated from new buffer function.

Square, Octagon, etc. generated from new buffer function.

Zoom in of different buffer smoothnesses.

Zoom in of different buffer smoothnesses.

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