· 9 years ago · Dec 08, 2016, 02:36 AM
1CREATE TABLE my_polygon (
2 my_polygon_id SERIAL PRIMARY KEY,
3 common_id INTEGER NOT NULL,
4 value1 NUMERIC NOT NULL,
5 value2 NUMERIC NOT NULL,
6 value3 NUMERIC NOT NULL,
7 geom GEOMETRY(Polygon) NOT NULL
8)
9;
10
11CREATE INDEX ON my_polygon (common_id);
12CREATE INDEX ON my_polygon USING GIST (common_id, geom);
13
14CREATE TABLE my_point (
15 my_point_id SERIAL PRIMARY KEY,
16 common_id INTEGER NOT NULL,
17 pointvalue NUMERIC NOT NULL,
18 geom GEOMETRY(Point) NOT NULL
19);
20
21CREATE INDEX ON my_point (common_id);
22CREATE INDEX ON my_point USING GIST (common_id, geom);
23
24SELECT DISTINCT ON (my_point.my_point_id)
25 my_polygon.*,
26 my_point.my_point_id,
27 my_point.pointvalue,
28 my_point.geom AS pointgeom
29FROM my_polygon
30JOIN my_point ON my_point.common_id = my_polygon.common_id AND ST_Contains(my_polygon.geom, my_point.geom)
31WHERE my_polygon.common_id = 1
32ORDER BY my_point.my_point_id, my_polygon.my_polygon_id
33
34SELECT *
35FROM (
36 SELECT DISTINCT ON (my_point.my_point_id)
37 my_polygon.*,
38 my_point.my_point_id,
39 my_point.pointvalue,
40 my_point.geom AS pointgeom
41 FROM my_polygon
42 JOIN my_point ON my_point.common_id = my_polygon.common_id AND ST_Contains(my_polygon.geom, my_point.geom)
43 ORDER BY my_point.my_point_id, my_polygon.my_polygon_id
44) point_with_polygon
45
46CREATE EXTENSION IF NOT EXISTS postgis;
47CREATE EXTENSION IF NOT EXISTS btree_gist;
48
49
50-- DROP FUNCTION ST_GeneratePoints(geometry, numeric);
51DO $doblock$
52BEGIN
53 IF NOT EXISTS(SELECT * FROM pg_proc WHERE UPPER(proname) = UPPER('ST_GeneratePoints')) THEN
54 -- Create naive ST_GeneratePoints if version of PostGIS is not new enough
55 CREATE FUNCTION ST_GeneratePoints(g geometry, npoints numeric)
56 RETURNS geometry
57 VOLATILE
58 RETURNS NULL ON NULL INPUT
59 LANGUAGE plpgsql
60 AS $$
61 DECLARE
62 num_to_generate INTEGER := npoints::INTEGER;
63 x_min FLOAT := ST_XMin(g) + 0.0001;
64 x_max FLOAT := ST_XMax(g) - 0.0001;
65 y_min FLOAT := ST_YMin(g) + 0.0001;
66 y_max FLOAT := ST_YMax(g) - 0.0001;
67 temp_result GEOMETRY[];
68 result_array GEOMETRY[] := ARRAY[]::GEOMETRY[];
69 BEGIN
70 -- Reduce number of loops to reduce slow array_cat calls
71 WHILE num_to_generate > 0 LOOP
72 SELECT ARRAY_AGG(contained.point) INTO temp_result
73 FROM (
74 SELECT point
75 FROM (
76 SELECT ST_MakePoint(
77 x_min + random() * (x_max - x_min),
78 y_min + random() * (y_max - y_min)
79 ) point
80 -- Generate extras to reduce number of loops, at least 20
81 FROM generate_series(1, GREATEST(20, CEIL(1.5 * num_to_generate)))
82 ) candidate
83 WHERE ST_Contains(g, candidate.point)
84 -- Filter out extras if we have too many matches
85 LIMIT num_to_generate
86 ) contained
87 ;
88 IF ARRAY_LENGTH(temp_result, 1) > 0 THEN
89 result_array := array_cat(result_array, temp_result);
90 num_to_generate := npoints - COALESCE(ARRAY_LENGTH(result_array, 1), 0);
91 END IF;
92 END LOOP;
93 RETURN (SELECT ST_Union(point) FROM UNNEST(result_array) result (point));
94 END;
95 $$;
96 RAISE NOTICE 'Created ST_GeneratePoints';
97 ELSE
98 RAISE NOTICE 'ST_GeneratePoints exists';
99 END IF;
100END
101$doblock$
102;
103
104DROP TABLE IF EXISTS my_polygon;
105
106CREATE TABLE my_polygon (
107 my_polygon_id SERIAL PRIMARY KEY,
108 common_id INTEGER NOT NULL,
109 value1 NUMERIC NOT NULL,
110 value2 NUMERIC NOT NULL,
111 value3 NUMERIC NOT NULL,
112 geom GEOMETRY(Polygon) NOT NULL
113)
114;
115
116CREATE INDEX ON my_polygon (common_id);
117CREATE INDEX ON my_polygon USING GIST (common_id, geom);
118
119
120WITH common AS (
121 SELECT
122 common_id,
123 random() * 5000 AS common_x_translate,
124 random() * 5000 AS common_y_translate
125 FROM (
126 SELECT TRUNC(random() * 1000) + 1 AS common_id
127 FROM generate_series(1, 100)
128 UNION
129 SELECT 1
130 ) a
131),
132geom_set_with_small_overlaps AS (
133 SELECT
134 ST_MakeEnvelope(
135 x.translate,
136 y.translate,
137 x.translate + 1.1,
138 y.translate + 1.1
139 ) AS geom
140 FROM
141 generate_series(0, 9) x (translate),
142 generate_series(0, 9) y (translate)
143)
144INSERT INTO my_polygon (common_id, value1, value2, value3, geom)
145SELECT
146 common_id,
147 random() * 100,
148 random() * 100,
149 random() * 100,
150 ST_Translate(geom, common_x_translate, common_y_translate)
151FROM common, geom_set_with_small_overlaps
152;
153
154DROP TABLE IF EXISTS my_point;
155
156CREATE TABLE my_point (
157 my_point_id SERIAL PRIMARY KEY,
158 common_id INTEGER NOT NULL,
159 pointvalue NUMERIC NOT NULL,
160 geom GEOMETRY(Point) NOT NULL
161);
162
163INSERT INTO my_point (common_id, pointvalue, geom)
164SELECT
165 common_id,
166 random() * 100,
167 (ST_Dump(ST_GeneratePoints(extent, FLOOR(5000 + random() * 15000)::NUMERIC))).geom
168FROM (
169 SELECT
170 common_id,
171 -- Small negative buffer prevents lying on the outer edge
172 ST_Buffer(ST_Extent(geom), - 0.0001) AS extent
173 FROM my_polygon
174 GROUP BY common_id
175) common
176UNION ALL
177SELECT
178 common_id,
179 random() * 100,
180 (ST_Dump(ST_GeneratePoints(intersection, TRUNC(random() * 5)::NUMERIC))).geom
181FROM (
182 SELECT
183 p1.common_id,
184 p1.my_polygon_id AS id1,
185 p2.my_polygon_id AS id2,
186 ST_Intersection(p1.geom, p2.geom) AS intersection
187 FROM my_polygon p1
188 JOIN my_polygon p2 ON (
189 p1.my_polygon_id < p2.my_polygon_id AND
190 p1.common_id = p2.common_id AND
191 ST_Intersects(p1.geom, p2.geom)
192 )
193) a
194;
195
196CREATE INDEX ON my_point (common_id);
197CREATE INDEX ON my_point USING GIST (common_id, geom);