postgis2.2_9.5_multiple_indexes_problem.sql

DROP TABLE IF EXISTS regions_medial_subset CASCADE;
CREATE TABLE regions_medial_subset AS 
SELECT gid,
	ST_ApproximateMedialAxis(ST_Buffer(geom, 0)) AS geom,
	scalerank, featurecla, "name", namealt, region, subregion	
	FROM regions_polys_subset
		WHERE
	"name" != 'AUSTRALIA' AND "name" != 'NORTH AMERICA' AND
	"name" != 'SOUTH AMERICA' AND "name" != 'EUROPE' AND
	"name" !=  'ASIA' AND "name" != 'AFRICA' AND "name" != 'ANTARCTICA'
	;

--SELECT GeometryType(geom) FROM regions_medial LIMIT 1;

DROP TABLE IF EXISTS regions_points_subset CASCADE;
CREATE TABLE regions_points_subset AS
SELECT gid AS oid,
	(ST_DumpPoints(ST_Segmentize(geom,1000))).geom AS geom,
	scalerank, featurecla, "name", namealt, region, subregion
	FROM regions_medial_subset;

ALTER TABLE regions_points_subset ADD COLUMN gid serial NOT NULL;
ALTER TABLE regions_points_subset ADD PRIMARY KEY (gid);
CREATE INDEX regions_points_subset_geom_idx
  ON regions_points_subset
  USING GIST (geom); 

DROP TABLE IF EXISTS regions_dpoints_subset CASCADE;
CREATE TABLE regions_dpoints_subset AS
SELECT gid,
	(ST_DumpPoints(ST_Segmentize(geom,500))).geom AS geom,
	scalerank, featurecla, "name", namealt, region, subregion
	FROM regions_polys_subset;

CREATE INDEX regions_dpoints_subset_geom_idx
  ON regions_dpoints_subset
  USING GIST (geom);	

CREATE OR REPLACE FUNCTION zz_1nn_d(geometry) RETURNS float AS $$

-- Here are my wonderful points to KNN search:
WITH index_query AS (

	SELECT edge.geom AS geom
	FROM (SELECT * FROM regions_dpoints WHERE gid = 881) AS edge
-- This is my query point
	ORDER BY $1
		<->
	edge.geom LIMIT 1
	)
SELECT ST_Distance($1,geom) FROM index_query

$$ LANGUAGE SQL;

CREATE OR REPLACE FUNCTION zz_2nn_angle(geometry) RETURNS float AS $$

-- Here are my wonderful points to KNN search:
WITH index_query AS (

	SELECT edge.geom AS geom
	FROM (SELECT * FROM regions_dpoints_subset) AS edge
-- This is my query point
	ORDER BY $1
		<->
	edge.geom LIMIT 2
	),
templine AS (
	SELECT ST_MakeLine(geom) AS geom FROM index_query
),
angle1 AS (
	SELECT ST_Azimuth(ST_StartPoint(geom), $1) angle FROM templine
),
angle2 AS (
	SELECT ST_Azimuth(ST_EndPoint(geom), $1) angle FROM templine
)
SELECT a1.angle - a2.angle FROM angle1 a1, angle2 a2

$$ LANGUAGE SQL;

DROP TABLE IF EXISTS test_angle_subset CASCADE;
CREATE TABLE test_angle_subset AS
WITH returnline AS (
	SELECT oid, gid, geom,
		zz_1nn_d(geom) AS distance,
		abs(degrees(zz_2nn_angle(geom))) AS angle FROM
		(SELECT * FROM regions_points_subset) subset
)
SELECT * FROM returnline
;

SELECT postgis_full_version();
SELECT ST_AsEWKT(geom) FROM regions_dpoints LIMIT 1;

ALTER EXTENSION postgis UPDATE TO "2.2.1dev";