pgLatLon
changeset 0:3b70e93cc07d v0.1
Version 0.1 (initial commit)
| author | jbe |
|---|---|
| date | Sun Aug 21 17:43:48 2016 +0200 (2016-08-21) |
| parents | |
| children | 2393084ef356 |
| files | GNUmakefile LICENSE Makefile README.html README.mkd create_test_db.data.sql create_test_db.schema.sql create_test_db.sh latlon--0.1.sql latlon-v0001.c latlon.control make-doc.sh |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/GNUmakefile Sun Aug 21 17:43:48 2016 +0200 1.3 @@ -0,0 +1,7 @@ 1.4 +EXTENSION = latlon 1.5 +DATA = latlon--0.1.sql 1.6 +MODULES = latlon-v0001 1.7 + 1.8 +PG_CONFIG = pg_config 1.9 +PGXS := $(shell $(PG_CONFIG) --pgxs) 1.10 +include $(PGXS)
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 2.2 +++ b/LICENSE Sun Aug 21 17:43:48 2016 +0200 2.3 @@ -0,0 +1,19 @@ 2.4 +Copyright (c) 2016 Public Software Group e. V., Berlin, Germany 2.5 + 2.6 +Permission is hereby granted, free of charge, to any person obtaining a 2.7 +copy of this software and associated documentation files (the "Software"), 2.8 +to deal in the Software without restriction, including without limitation 2.9 +the rights to use, copy, modify, merge, publish, distribute, sublicense, 2.10 +and/or sell copies of the Software, and to permit persons to whom the 2.11 +Software is furnished to do so, subject to the following conditions: 2.12 + 2.13 +The above copyright notice and this permission notice shall be included in 2.14 +all copies or substantial portions of the Software. 2.15 + 2.16 +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 2.17 +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 2.18 +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 2.19 +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 2.20 +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 2.21 +FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 2.22 +DEALINGS IN THE SOFTWARE.
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 3.2 +++ b/Makefile Sun Aug 21 17:43:48 2016 +0200 3.3 @@ -0,0 +1,8 @@ 3.4 +all:: 3.5 + gmake all 3.6 +clean:: 3.7 + gmake clean 3.8 +install:: 3.9 + gmake install 3.10 +uninstall:: 3.11 + gmake uninstall
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 4.2 +++ b/README.html Sun Aug 21 17:43:48 2016 +0200 4.3 @@ -0,0 +1,450 @@ 4.4 +<html><head><title>pgLatLon v0.1 documentation</title></head><body> 4.5 +<h1>pgLatLon v0.1 documentation</h1> 4.6 + 4.7 +<p>pgLatLon is a spatial database extension for the PostgreSQL object-relational 4.8 +database management system providing geographic data types and spatial indexing 4.9 +for the WGS-84 spheroid.</p> 4.10 + 4.11 +<p>While many other spatial databases still use imprecise bounding boxes for many 4.12 +operations, pgLatLon supports more precise geometric calculations for all 4.13 +implemented operators. Efficient indexing of geometric objects is provided 4.14 +using fractal indices. Optimizations on bit level (including logarithmic 4.15 +compression) allow for a highly memory-efficient non-overlapping index suitable 4.16 +for huge datasets.</p> 4.17 + 4.18 +<p>Unlike competing spatial extensions for PostgreSQL, pgLatLon is available under 4.19 +the permissive MIT/X11 license to avoid problems with viral licenses like the 4.20 +GPLv2/v3.</p> 4.21 + 4.22 +<h2>Installation</h2> 4.23 + 4.24 +<h3>Automatic installation</h3> 4.25 + 4.26 +<p>Prerequisites:</p> 4.27 + 4.28 +<ul> 4.29 +<li>Ensure that the <code>pg_config</code> binary is in your path (shipped with PostgreSQL).</li> 4.30 +<li>Ensure that GNU Make is available (either as <code>make</code> or <code>gmake</code>).</li> 4.31 +</ul> 4.32 + 4.33 +<p>Then simply type:</p> 4.34 + 4.35 +<pre><code>make install 4.36 +</code></pre> 4.37 + 4.38 +<h3>Manual installation</h3> 4.39 + 4.40 +<p>It is also possible to compile and install the extension without GNU Make as 4.41 +follows:</p> 4.42 + 4.43 +<pre><code>cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0001.so latlon-v0001.c 4.44 +cp latlon-v0001.so `pg_config --pkglibdir` 4.45 +cp latlon.control `pg_config --sharedir`/extension/ 4.46 +cp latlon--0.1.sql `pg_config --sharedir`/extension/ 4.47 +</code></pre> 4.48 + 4.49 +<h3>Loading the extension</h3> 4.50 + 4.51 +<p>After installation, you can create a database and load the extension as 4.52 +follows:</p> 4.53 + 4.54 +<pre><code>% createdb test_database 4.55 +% psql test_database 4.56 +psql (9.5.4) 4.57 +Type "help" for help. 4.58 + 4.59 +test_database=# CREATE EXTENSION latlon; 4.60 +</code></pre> 4.61 + 4.62 +<h2>Reference</h2> 4.63 + 4.64 +<h3>1. Types</h3> 4.65 + 4.66 +<p>pgLatLon provides four geographic types: <code>epoint</code>, <code>ebox</code>, <code>ecircle</code>, and 4.67 +<code>ecluster</code>.</p> 4.68 + 4.69 +<h4><code>epoint</code></h4> 4.70 + 4.71 +<p>A point on the earth spheroid (WGS-84).</p> 4.72 + 4.73 +<p>The text input format is <code>'[N|S]<float> [E|W]<float>'</code>, where each float is in 4.74 +degrees. Note the required white space between the latitude and longitude 4.75 +components. Each floating point number may have a sign, in which case <code>N</code>/<code>S</code> 4.76 +or <code>E</code>/<code>W</code> are switched respectively (e.g. <code>E-5</code> is the same as <code>W5</code>).</p> 4.77 + 4.78 +<p>An <code>epoint</code> may also be created from two floating point numbers by calling 4.79 +<code>epoint(latitude, longitude)</code>, where positive latitudes are used for the 4.80 +northern hemisphere, negative latitudes are used for the southern hemisphere, 4.81 +positive longitudes indicate positions east of the prime meridian, and negative 4.82 +longitudes indicate positions west of the prime meridian.</p> 4.83 + 4.84 +<p>Latitudes exceeding -90 or +90 degrees are truncated to -90 or +90 4.85 +respectively, in which case a warning will be issued. Longitudes exceeding -180 4.86 +or +180 degrees will be converted to values between -180 and +180 (both 4.87 +inclusive) by adding or substracting a multiple of 360 degrees, in which case a 4.88 +notice will be issued.</p> 4.89 + 4.90 +<p>If the latitude is -90 or +90 (south pole or north pole), a longitude value is 4.91 +still stored in the datum, and if a point is on the prime meridian or the 4.92 +180th meridian, the east/west bit is also stored in the datum. In case of the 4.93 +prime meridian, this is done by storing a floating point value of -0 for 4.94 +0 degrees west and a value of +0 for 0 degrees east. In case of the 4.95 +180th meridian, this is done by storing -180 or +180 respectively. The equality 4.96 +operator, however, returns true when the same points on earth are described, 4.97 +i.e. the longitude is ignored for the poles, and 180 degrees west is considered 4.98 +to be equal to 180 degrees east.</p> 4.99 + 4.100 +<h4><code>ebox</code></h4> 4.101 + 4.102 +<p>An area on earth demarcated by a southern and northern latitude, and a western 4.103 +and eastern longitude (all given in WGS-84).</p> 4.104 + 4.105 +<p>The text input format is 4.106 +<code>'{N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float>'</code>, where each float is in 4.107 +degrees. The ordering of the four white-space separated blocks is not 4.108 +significant. To include the 180th meridian, one longitude boundary must be 4.109 +equal to or exceed <code>W180</code> or <code>E180</code>, e.g. <code>'N10 N20 E170 E190'</code>.</p> 4.110 + 4.111 +<p>A special value is the empty area, denoted by the text represenation <code>'empty'</code>. 4.112 +Such an <code>ebox</code> does not contain any point.</p> 4.113 + 4.114 +<p>An <code>ebox</code> may also be created from four floating point numbers by calling 4.115 +<code>ebox(min_latitude, max_latitude, min_longitude, max_longitude)</code>, where 4.116 +positive values are used for north and east, and negative values are used for 4.117 +south and west. If <code>min_latitude</code> is strictly greater than <code>max_latitude</code>, an 4.118 +empty <code>ebox</code> is created. If <code>min_longitude</code> is greater than <code>max_longitude</code> and 4.119 +if both longitudes are between -180 and +180 degrees, then the area oriented in 4.120 +such way that the 180th meridian is included.</p> 4.121 + 4.122 +<p>If the longitude span is less than 120 degrees, an <code>ebox</code> may be alternatively 4.123 +created from two <code>epoints</code> in the following way: <code>ebox(epoint(lat1, lon1), 4.124 +epoint(lat2, lon2))</code>. In this case <code>lat1</code> and <code>lat2</code> as well as <code>lon1</code> and 4.125 +<code>lon2</code> can be swapped without any impact.</p> 4.126 + 4.127 +<h4><code>ecircle</code></h4> 4.128 + 4.129 +<p>An area containing all points not farther away from a given center point 4.130 +(WGS-84) than a given radius.</p> 4.131 + 4.132 +<p>The text input format is <code>'{N|S}<float> {E|W}<float> <float>'</code>, where the first 4.133 +two floats denote the center point in degrees and the third float denotes the 4.134 +radius in meters. A radius equal to minus infinity denotes an empty circle 4.135 +which contains no point at all (despite having a center), while a radius equal 4.136 +to zero denotes a circle that includes a single point.</p> 4.137 + 4.138 +<p>An <code>ecircle</code> may also be created by calling <code>ecircle(epoint(...), radius)</code> or 4.139 +from three floating point numbers by calling <code>ecircle(latitude, longitude, 4.140 +radius)</code>.</p> 4.141 + 4.142 +<h4><code>ecluster</code></h4> 4.143 + 4.144 +<p>A collection of points, paths, polygons, and outlines on the WGS-84 spheroid. 4.145 +Each path, polygon, or outline must cover a longitude range of less than 4.146 +180 degrees to avoid ambiguities.</p> 4.147 + 4.148 +<p>The text input format is a white-space separated list of the following items:</p> 4.149 + 4.150 +<ul> 4.151 +<li><code>point ({N|S}<float> {E|W}<float>)</code></li> 4.152 +<li><code>path ({N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> ...)</code></li> 4.153 +<li><code>outline ({N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> ...)</code></li> 4.154 +<li><code>polygon ({N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> ...)</code></li> 4.155 +</ul> 4.156 + 4.157 +<p>Paths are open by default (i.e. there is no connection from the last point in 4.158 +the list to the first point in the list). Outlines and polygons, in contrast, 4.159 +are automatically closed (i.e. there is a line segment from the last point in 4.160 +the list to the first point in the list) which means the first point should not 4.161 +be repeated as last point in the list. Polygons are filled, outlines are not.</p> 4.162 + 4.163 +<h3>2. Indices</h3> 4.164 + 4.165 +<p>Two kinds of indices are supported: B-tree and GiST indices.</p> 4.166 + 4.167 +<h4>B-tree indices</h4> 4.168 + 4.169 +<p>A B-tree index can be used for simple equality searches and is supported by the 4.170 +<code>epoint</code>, <code>ebox</code>, and <code>ecircle</code> data types. B-tree indices can not be used for 4.171 +geographic searches.</p> 4.172 + 4.173 +<h4>GiST indices</h4> 4.174 + 4.175 +<p>For geographic searches, GiST indices must be used. The <code>epoint</code>, <code>ecircle</code>, 4.176 +and <code>ecluster</code> data types support GiST indexing. A GiST index for geographic 4.177 +searches can be created as follows:</p> 4.178 + 4.179 +<pre><code>CREATE TABLE tbl ( 4.180 + id serial4 PRIMARY KEY, 4.181 + loc epoint NOT NULL ); 4.182 + 4.183 +CREATE INDEX name_of_index ON tbl USING gist (loc); 4.184 +</code></pre> 4.185 + 4.186 +<p>GiST indices also support nearest neighbor searches when using the distance 4.187 +operator (<code><-></code>) in the ORDER BY clause.</p> 4.188 + 4.189 +<h4>Indices on other data types (e.g. GeoJSON)</h4> 4.190 + 4.191 +<p>Note that further types can be indexed by using an index on an expression with 4.192 +a conversion function. One conversion function provided by pgLatLon is the 4.193 +<code>GeoJSON_to_ecluster(float8, float8, text)</code> function:</p> 4.194 + 4.195 +<pre><code>CREATE TABLE tbl ( 4.196 + id serial4 PRIMARY KEY, 4.197 + loc jsonb NOT NULL ); 4.198 + 4.199 +CREATE INDEX name_of_index ON tbl USING gist((GeoJSON_to_ecluster("loc"))); 4.200 +</code></pre> 4.201 + 4.202 +<p>When using the conversion function in an expression, the index will be used 4.203 +automatically:</p> 4.204 + 4.205 +<pre><code>SELECT * FROM tbl WHERE GeoJSON_to_ecluster("loc") && 'N50 E10 10000'::ecircle; 4.206 +</code></pre> 4.207 + 4.208 +<h3>3. Operators</h3> 4.209 + 4.210 +<h4>Equality operator <code>=</code></h4> 4.211 + 4.212 +<p>Tests if two geographic objects are equal.</p> 4.213 + 4.214 +<p>The longitude is ignored for the poles, and 180 degrees west is considered to 4.215 +be equal to 180 degrees east.</p> 4.216 + 4.217 +<p>For boxes and circles, two empty objects are considered equal. (Note that a 4.218 +circle is not empty if the radius is zero but only if it is negative infinity, 4.219 +i.e. smaller than zero.) Two circles with a positive infinite radius are also 4.220 +considered equal.</p> 4.221 + 4.222 +<p>Implemented for:</p> 4.223 + 4.224 +<ul> 4.225 +<li><code>epoint = epoint</code></li> 4.226 +<li><code>ebox = ebox</code></li> 4.227 +<li><code>ecircle = ecircle</code></li> 4.228 +</ul> 4.229 + 4.230 +<p>The negation is the inequality operator (<code><></code> or <code>!=</code>).</p> 4.231 + 4.232 +<h4>Linear ordering operators <code><<<</code>, <code><<<=</code>, <code>>>>=</code>, <code>>>></code></h4> 4.233 + 4.234 +<p>These operators create an arbitrary (but well-defined) linear ordering of 4.235 +geographic objects, which is used internally for B-tree indexing and merge 4.236 +joins. These operators will usually not be used by an application programmer.</p> 4.237 + 4.238 +<h4>Overlap operator <code>&&</code></h4> 4.239 + 4.240 +<p>Tests if two geographic objects have at least one point in common. Currently 4.241 +implemented for:</p> 4.242 + 4.243 +<ul> 4.244 +<li><code>epoint && ebox</code></li> 4.245 +<li><code>epoint && ecircle</code></li> 4.246 +<li><code>epoint && ecluster</code></li> 4.247 +<li><code>ebox && ebox</code></li> 4.248 +<li><code>ecircle && ecircle</code></li> 4.249 +<li><code>ecircle && ecluster</code></li> 4.250 +</ul> 4.251 + 4.252 +<p>The <code>&&</code> operator is commutative, i.e. <code>a && b</code> is the same as <code>b && a</code>. Each 4.253 +commutation is supported as well.</p> 4.254 + 4.255 +<h4>Distance operator <code><-></code></h4> 4.256 + 4.257 +<p>Calculates the shortest distance between two geographic objects in meters (zero 4.258 +if the objects are overlapping). Currently implemented for:</p> 4.259 + 4.260 +<ul> 4.261 +<li><code>epoint <-> epoint</code></li> 4.262 +<li><code>epoint <-> ecircle</code></li> 4.263 +<li><code>epoint <-> ecluster</code></li> 4.264 +<li><code>ecircle <-> ecircle</code></li> 4.265 +<li><code>ecircle <-> ecluster</code></li> 4.266 +</ul> 4.267 + 4.268 +<p>The <code><-></code> operator is commutative, i.e. <code>a <-> b</code> is the same as <code>b <-> a</code>. 4.269 +Each commutation is supported as well.</p> 4.270 + 4.271 +<p>For short distances, the result is very accurate (i.e. respects the dimensions 4.272 +of the WGS-84 spheroid). For longer distances in the order of magnitude of 4.273 +earth's radius or greater, the value is only approximate (but the error is 4.274 +still less than 0.2% as long as no polygons with very long edges are involved).</p> 4.275 + 4.276 +<p>The functions <code>distance(epoint, epoint)</code> and <code>distance(ecluster, epoint)</code> can 4.277 +be used as an alias for this operator.</p> 4.278 + 4.279 +<p>Note: In case of radial searches with a fixed radius, this operator should 4.280 +not be used. Instead, an <code>ecircle</code> should be created and used in combination 4.281 +with the overlap operator (<code>&&</code>). Alternatively, the functions 4.282 +<code>distance_within(epoint, epoint, float8)</code> or <code>distance_within(ecluster, epoint, 4.283 +float8)</code> can be used for fixed-radius searches.</p> 4.284 + 4.285 +<h3>4. Functions</h3> 4.286 + 4.287 +<h4><code>center(circle)</code></h4> 4.288 + 4.289 +<p>Returns the center of an <code>ecircle</code> as an <code>epoint</code>.</p> 4.290 + 4.291 +<h4><code>distance(epoint, epoint)</code></h4> 4.292 + 4.293 +<p>Calculates the distance between two <code>epoint</code> datums in meters. This function is 4.294 +an alias for the distance operator <code><-></code>.</p> 4.295 + 4.296 +<p>Note: In case of radial searches with a fixed radius, this function should not be 4.297 +used. Use <code>distance_within(epoint, epoint, float8)</code> instead.</p> 4.298 + 4.299 +<h4><code>distance(ecluster, epoint)</code></h4> 4.300 + 4.301 +<p>Calculates the distance from an <code>ecluster</code> to an <code>epoint</code> in meters. This 4.302 +function is an alias for the distance operator <code><-></code>.</p> 4.303 + 4.304 +<p>Note: In case of radial searches with a fixed radius, this function should not be 4.305 +used. Use <code>distance_within(epoint, epoint, float8)</code> instead.</p> 4.306 + 4.307 +<h4><code>distance_within(</code>variable <code>epoint,</code> fixed <code>epoint,</code> radius <code>float8)</code></h4> 4.308 + 4.309 +<p>Checks if the distance between two <code>epoint</code> datums is not greater than a given 4.310 +value (search radius).</p> 4.311 + 4.312 +<p>Note: In case of radial searches with a fixed radius, the first argument must 4.313 +be used for the table column, while the second argument must be used for the 4.314 +search center. Otherwise an existing index cannot be used.</p> 4.315 + 4.316 +<h4><code>distance_within(</code>variable <code>ecluster,</code> fixed <code>epoint,</code> radius <code>float8)</code></h4> 4.317 + 4.318 +<p>Checks if the distance from an <code>ecluster</code> to an <code>epoint</code> is not greater than a 4.319 +given value (search radius).</p> 4.320 + 4.321 +<h4><code>ebox(</code>latmin <code>float8,</code> latmax <code>float8,</code> lonmin <code>float8,</code> lonmax <code>float8)</code></h4> 4.322 + 4.323 +<p>Creates a new <code>ebox</code> with the given boundaries. 4.324 +See "1. Types", subsection <code>ebox</code> for details.</p> 4.325 + 4.326 +<h4><code>ebox(epoint, epoint)</code></h4> 4.327 + 4.328 +<p>Creates a new <code>ebox</code>. This function may only be used if the longitude 4.329 +difference is less than or equal to 120 degrees. 4.330 +See "1. Types", subsection <code>ebox</code> for details.</p> 4.331 + 4.332 +<h4><code>ecircle(epoint, float8)</code></h4> 4.333 + 4.334 +<p>Creates an <code>ecircle</code> with the given center point and radius.</p> 4.335 + 4.336 +<h4><code>ecircle(</code>latitude <code>float8,</code> longitude <code>float8,</code> radius <code>float8)</code></h4> 4.337 + 4.338 +<p>Creates an <code>ecircle</code> with the given center point and radius.</p> 4.339 + 4.340 +<h4><code>ecluster_concat(ecluster, ecluster)</code></h4> 4.341 + 4.342 +<p>Combines two clusters to form a new <code>ecluster</code> by uniting all entries of both 4.343 +clusters. Note that two overlapping areas of polygons annihilate each other 4.344 +(which may be used to create polygons with holes).</p> 4.345 + 4.346 +<h4><code>ecluster_concat(ecluster[])</code></h4> 4.347 + 4.348 +<p>Creates a new <code>ecluster</code> that unites all entries of all clusters in the passed 4.349 +array. Note that two overlapping areas of polygons annihilate each other (which 4.350 +may be used to create polygons with holes).</p> 4.351 + 4.352 +<h4><code>ecluster_create_multipoint(epoint[])</code></h4> 4.353 + 4.354 +<p>Creates a new <code>ecluster</code> which contains multiple points.</p> 4.355 + 4.356 +<h4><code>ecluster_create_outline(epoint[])</code></h4> 4.357 + 4.358 +<p>Creates a new <code>ecluster</code> that is an outline given by the passed points.</p> 4.359 + 4.360 +<h4><code>ecluster_create_path(epoint[])</code></h4> 4.361 + 4.362 +<p>Creates a new <code>ecluster</code> that is a path given by the passed points.</p> 4.363 + 4.364 +<h4><code>ecluster_create_polygon(epoint[])</code></h4> 4.365 + 4.366 +<p>Creates a new <code>ecluster</code> that is a polygon given by the passed points.</p> 4.367 + 4.368 +<h4><code>ecluster_extract_outlines(ecluster)</code></h4> 4.369 + 4.370 +<p>Set-returning function that returns the outlines of an <code>ecluster</code> as <code>epoint[]</code> 4.371 +rows.</p> 4.372 + 4.373 +<h4><code>ecluster_extract_paths(ecluster)</code></h4> 4.374 + 4.375 +<p>Set-returning function that returns the paths of an <code>ecluster</code> as <code>epoint[]</code> 4.376 +rows.</p> 4.377 + 4.378 +<h4><code>ecluster_extract_points(ecluster)</code></h4> 4.379 + 4.380 +<p>Set-returning function that returns the points of an <code>ecluster</code> as <code>epoint</code> 4.381 +rows.</p> 4.382 + 4.383 +<h4><code>ecluster_extract_polygons(ecluster)</code></h4> 4.384 + 4.385 +<p>Set-returning function that returns the polygons of an <code>ecluster</code> as <code>epoint[]</code> 4.386 +rows.</p> 4.387 + 4.388 +<h4><code>empty_ebox</code>()</h4> 4.389 + 4.390 +<p>Returns the empty <code>ebox</code>. 4.391 +See "1. Types", subsection <code>ebox</code> for details.</p> 4.392 + 4.393 +<h4><code>epoint(</code>latitude <code>float8,</code> longitude <code>float8)</code></h4> 4.394 + 4.395 +<p>Returns an <code>epoint</code> with the given latitude and longitude.</p> 4.396 + 4.397 +<h4><code>epoint_latlon(</code>latitude <code>float8,</code> longitude <code>float8)</code></h4> 4.398 + 4.399 +<p>Alias for <code>epoint(float8, float8)</code>.</p> 4.400 + 4.401 +<h4><code>epoint_lonlat(</code>longitude <code>float8,</code> latitude <code>float8)</code></h4> 4.402 + 4.403 +<p>Same as <code>epoint(float8, float8)</code> but with arguments reversed.</p> 4.404 + 4.405 +<h4><code>GeoJSON_to_epoint(jsonb, text)</code></h4> 4.406 + 4.407 +<p>Maps a GeoJSON object of type "Point" or "Feature" (which contains a 4.408 +"Point") to an <code>epoint</code> datum. For any other JSON objects, NULL is returned.</p> 4.409 + 4.410 +<p>The second parameter (which defaults to <code>epoint_lonlat</code>) may be set to a name 4.411 +of a conversion function that transforms two coordinates (two <code>float8</code> 4.412 +parameters) to an <code>epoint</code>.</p> 4.413 + 4.414 +<h4><code>GeoJSON_to_ecluster(jsonb, text)</code></h4> 4.415 + 4.416 +<p>Maps a (valid) GeoJSON object to an <code>ecluster</code>. Note that this function 4.417 +does not check whether the JSONB object is a valid GeoJSON object.</p> 4.418 + 4.419 +<p>The second parameter (which defaults to <code>epoint_lonlat</code>) may be set to a name 4.420 +of a conversion function that transforms two coordinates (two <code>float8</code> 4.421 +parameters) to an <code>epoint</code>.</p> 4.422 + 4.423 +<h4><code>max_latitude(ebox)</code></h4> 4.424 + 4.425 +<p>Returns the northern boundary of a given <code>ebox</code> in degrees between -90 and +90.</p> 4.426 + 4.427 +<h4><code>max_longitude(ebox)</code></h4> 4.428 + 4.429 +<p>Returns the eastern boundary of a given <code>ebox</code> in degrees between -180 and +180 4.430 +(both inclusive).</p> 4.431 + 4.432 +<h4><code>min_latitude(ebox)</code></h4> 4.433 + 4.434 +<p>Returns the southern boundary of a given <code>ebox</code> in degrees between -90 and +90.</p> 4.435 + 4.436 +<h4><code>min_longitude(ebox)</code></h4> 4.437 + 4.438 +<p>Returns the western boundary of a given <code>ebox</code> in degrees between -180 and +180 4.439 +(both inclusive).</p> 4.440 + 4.441 +<h4><code>latitude(epoint)</code></h4> 4.442 + 4.443 +<p>Returns the latitude value of an <code>epoint</code> in degrees between -90 and +90.</p> 4.444 + 4.445 +<h4><code>longitude(epoint)</code></h4> 4.446 + 4.447 +<p>Returns the longitude value of an <code>epoint</code> in degrees between -180 and +180 4.448 +(both inclusive).</p> 4.449 + 4.450 +<h4><code>radius(ecircle)</code></h4> 4.451 + 4.452 +<p>Returns the radius of an <code>ecircle</code> in meters.</p> 4.453 +</body></html>
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 5.2 +++ b/README.mkd Sun Aug 21 17:43:48 2016 +0200 5.3 @@ -0,0 +1,438 @@ 5.4 +pgLatLon v0.1 documentation 5.5 +=========================== 5.6 + 5.7 +pgLatLon is a spatial database extension for the PostgreSQL object-relational 5.8 +database management system providing geographic data types and spatial indexing 5.9 +for the WGS-84 spheroid. 5.10 + 5.11 +While many other spatial databases still use imprecise bounding boxes for many 5.12 +operations, pgLatLon supports more precise geometric calculations for all 5.13 +implemented operators. Efficient indexing of geometric objects is provided 5.14 +using fractal indices. Optimizations on bit level (including logarithmic 5.15 +compression) allow for a highly memory-efficient non-overlapping index suitable 5.16 +for huge datasets. 5.17 + 5.18 +Unlike competing spatial extensions for PostgreSQL, pgLatLon is available under 5.19 +the permissive MIT/X11 license to avoid problems with viral licenses like the 5.20 +GPLv2/v3. 5.21 + 5.22 + 5.23 +Installation 5.24 +------------ 5.25 + 5.26 +### Automatic installation 5.27 + 5.28 +Prerequisites: 5.29 + 5.30 +* Ensure that the `pg_config` binary is in your path (shipped with PostgreSQL). 5.31 +* Ensure that GNU Make is available (either as `make` or `gmake`). 5.32 + 5.33 +Then simply type: 5.34 + 5.35 + make install 5.36 + 5.37 +### Manual installation 5.38 + 5.39 +It is also possible to compile and install the extension without GNU Make as 5.40 +follows: 5.41 + 5.42 + cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0001.so latlon-v0001.c 5.43 + cp latlon-v0001.so `pg_config --pkglibdir` 5.44 + cp latlon.control `pg_config --sharedir`/extension/ 5.45 + cp latlon--0.1.sql `pg_config --sharedir`/extension/ 5.46 + 5.47 +### Loading the extension 5.48 + 5.49 +After installation, you can create a database and load the extension as 5.50 +follows: 5.51 + 5.52 + % createdb test_database 5.53 + % psql test_database 5.54 + psql (9.5.4) 5.55 + Type "help" for help. 5.56 + 5.57 + test_database=# CREATE EXTENSION latlon; 5.58 + 5.59 + 5.60 +Reference 5.61 +--------- 5.62 + 5.63 +### 1. Types 5.64 + 5.65 +pgLatLon provides four geographic types: `epoint`, `ebox`, `ecircle`, and 5.66 +`ecluster`. 5.67 + 5.68 +#### `epoint` 5.69 + 5.70 +A point on the earth spheroid (WGS-84). 5.71 + 5.72 +The text input format is `'[N|S]<float> [E|W]<float>'`, where each float is in 5.73 +degrees. Note the required white space between the latitude and longitude 5.74 +components. Each floating point number may have a sign, in which case `N`/`S` 5.75 +or `E`/`W` are switched respectively (e.g. `E-5` is the same as `W5`). 5.76 + 5.77 +An `epoint` may also be created from two floating point numbers by calling 5.78 +`epoint(latitude, longitude)`, where positive latitudes are used for the 5.79 +northern hemisphere, negative latitudes are used for the southern hemisphere, 5.80 +positive longitudes indicate positions east of the prime meridian, and negative 5.81 +longitudes indicate positions west of the prime meridian. 5.82 + 5.83 +Latitudes exceeding -90 or +90 degrees are truncated to -90 or +90 5.84 +respectively, in which case a warning will be issued. Longitudes exceeding -180 5.85 +or +180 degrees will be converted to values between -180 and +180 (both 5.86 +inclusive) by adding or substracting a multiple of 360 degrees, in which case a 5.87 +notice will be issued. 5.88 + 5.89 +If the latitude is -90 or +90 (south pole or north pole), a longitude value is 5.90 +still stored in the datum, and if a point is on the prime meridian or the 5.91 +180th meridian, the east/west bit is also stored in the datum. In case of the 5.92 +prime meridian, this is done by storing a floating point value of -0 for 5.93 +0 degrees west and a value of +0 for 0 degrees east. In case of the 5.94 +180th meridian, this is done by storing -180 or +180 respectively. The equality 5.95 +operator, however, returns true when the same points on earth are described, 5.96 +i.e. the longitude is ignored for the poles, and 180 degrees west is considered 5.97 +to be equal to 180 degrees east. 5.98 + 5.99 +#### `ebox` 5.100 + 5.101 +An area on earth demarcated by a southern and northern latitude, and a western 5.102 +and eastern longitude (all given in WGS-84). 5.103 + 5.104 +The text input format is 5.105 +`'{N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float>'`, where each float is in 5.106 +degrees. The ordering of the four white-space separated blocks is not 5.107 +significant. To include the 180th meridian, one longitude boundary must be 5.108 +equal to or exceed `W180` or `E180`, e.g. `'N10 N20 E170 E190'`. 5.109 + 5.110 +A special value is the empty area, denoted by the text represenation `'empty'`. 5.111 +Such an `ebox` does not contain any point. 5.112 + 5.113 +An `ebox` may also be created from four floating point numbers by calling 5.114 +`ebox(min_latitude, max_latitude, min_longitude, max_longitude)`, where 5.115 +positive values are used for north and east, and negative values are used for 5.116 +south and west. If `min_latitude` is strictly greater than `max_latitude`, an 5.117 +empty `ebox` is created. If `min_longitude` is greater than `max_longitude` and 5.118 +if both longitudes are between -180 and +180 degrees, then the area oriented in 5.119 +such way that the 180th meridian is included. 5.120 + 5.121 +If the longitude span is less than 120 degrees, an `ebox` may be alternatively 5.122 +created from two `epoints` in the following way: `ebox(epoint(lat1, lon1), 5.123 +epoint(lat2, lon2))`. In this case `lat1` and `lat2` as well as `lon1` and 5.124 +`lon2` can be swapped without any impact. 5.125 + 5.126 +#### `ecircle` 5.127 + 5.128 +An area containing all points not farther away from a given center point 5.129 +(WGS-84) than a given radius. 5.130 + 5.131 +The text input format is `'{N|S}<float> {E|W}<float> <float>'`, where the first 5.132 +two floats denote the center point in degrees and the third float denotes the 5.133 +radius in meters. A radius equal to minus infinity denotes an empty circle 5.134 +which contains no point at all (despite having a center), while a radius equal 5.135 +to zero denotes a circle that includes a single point. 5.136 + 5.137 +An `ecircle` may also be created by calling `ecircle(epoint(...), radius)` or 5.138 +from three floating point numbers by calling `ecircle(latitude, longitude, 5.139 +radius)`. 5.140 + 5.141 +#### `ecluster` 5.142 + 5.143 +A collection of points, paths, polygons, and outlines on the WGS-84 spheroid. 5.144 +Each path, polygon, or outline must cover a longitude range of less than 5.145 +180 degrees to avoid ambiguities. 5.146 + 5.147 +The text input format is a white-space separated list of the following items: 5.148 + 5.149 +* `point ({N|S}<float> {E|W}<float>)` 5.150 +* `path ({N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> ...)` 5.151 +* `outline ({N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> ...)` 5.152 +* `polygon ({N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> {N|S}<float> {E|W}<float> ...)` 5.153 + 5.154 +Paths are open by default (i.e. there is no connection from the last point in 5.155 +the list to the first point in the list). Outlines and polygons, in contrast, 5.156 +are automatically closed (i.e. there is a line segment from the last point in 5.157 +the list to the first point in the list) which means the first point should not 5.158 +be repeated as last point in the list. Polygons are filled, outlines are not. 5.159 + 5.160 +### 2. Indices 5.161 + 5.162 +Two kinds of indices are supported: B-tree and GiST indices. 5.163 + 5.164 +#### B-tree indices 5.165 + 5.166 +A B-tree index can be used for simple equality searches and is supported by the 5.167 +`epoint`, `ebox`, and `ecircle` data types. B-tree indices can not be used for 5.168 +geographic searches. 5.169 + 5.170 +#### GiST indices 5.171 + 5.172 +For geographic searches, GiST indices must be used. The `epoint`, `ecircle`, 5.173 +and `ecluster` data types support GiST indexing. A GiST index for geographic 5.174 +searches can be created as follows: 5.175 + 5.176 + CREATE TABLE tbl ( 5.177 + id serial4 PRIMARY KEY, 5.178 + loc epoint NOT NULL ); 5.179 + 5.180 + CREATE INDEX name_of_index ON tbl USING gist (loc); 5.181 + 5.182 +GiST indices also support nearest neighbor searches when using the distance 5.183 +operator (`<->`) in the ORDER BY clause. 5.184 + 5.185 +#### Indices on other data types (e.g. GeoJSON) 5.186 + 5.187 +Note that further types can be indexed by using an index on an expression with 5.188 +a conversion function. One conversion function provided by pgLatLon is the 5.189 +`GeoJSON_to_ecluster(float8, float8, text)` function: 5.190 + 5.191 + CREATE TABLE tbl ( 5.192 + id serial4 PRIMARY KEY, 5.193 + loc jsonb NOT NULL ); 5.194 + 5.195 + CREATE INDEX name_of_index ON tbl USING gist((GeoJSON_to_ecluster("loc"))); 5.196 + 5.197 +When using the conversion function in an expression, the index will be used 5.198 +automatically: 5.199 + 5.200 + SELECT * FROM tbl WHERE GeoJSON_to_ecluster("loc") && 'N50 E10 10000'::ecircle; 5.201 + 5.202 +### 3. Operators 5.203 + 5.204 +#### Equality operator `=` 5.205 + 5.206 +Tests if two geographic objects are equal. 5.207 + 5.208 +The longitude is ignored for the poles, and 180 degrees west is considered to 5.209 +be equal to 180 degrees east. 5.210 + 5.211 +For boxes and circles, two empty objects are considered equal. (Note that a 5.212 +circle is not empty if the radius is zero but only if it is negative infinity, 5.213 +i.e. smaller than zero.) Two circles with a positive infinite radius are also 5.214 +considered equal. 5.215 + 5.216 +Implemented for: 5.217 + 5.218 +* `epoint = epoint` 5.219 +* `ebox = ebox` 5.220 +* `ecircle = ecircle` 5.221 + 5.222 +The negation is the inequality operator (`<>` or `!=`). 5.223 + 5.224 +#### Linear ordering operators `<<<`, `<<<=`, `>>>=`, `>>>` 5.225 + 5.226 +These operators create an arbitrary (but well-defined) linear ordering of 5.227 +geographic objects, which is used internally for B-tree indexing and merge 5.228 +joins. These operators will usually not be used by an application programmer. 5.229 + 5.230 +#### Overlap operator `&&` 5.231 + 5.232 +Tests if two geographic objects have at least one point in common. Currently 5.233 +implemented for: 5.234 + 5.235 +* `epoint && ebox` 5.236 +* `epoint && ecircle` 5.237 +* `epoint && ecluster` 5.238 +* `ebox && ebox` 5.239 +* `ecircle && ecircle` 5.240 +* `ecircle && ecluster` 5.241 + 5.242 +The `&&` operator is commutative, i.e. `a && b` is the same as `b && a`. Each 5.243 +commutation is supported as well. 5.244 + 5.245 +#### Distance operator `<->` 5.246 + 5.247 +Calculates the shortest distance between two geographic objects in meters (zero 5.248 +if the objects are overlapping). Currently implemented for: 5.249 + 5.250 +* `epoint <-> epoint` 5.251 +* `epoint <-> ecircle` 5.252 +* `epoint <-> ecluster` 5.253 +* `ecircle <-> ecircle` 5.254 +* `ecircle <-> ecluster` 5.255 + 5.256 +The `<->` operator is commutative, i.e. `a <-> b` is the same as `b <-> a`. 5.257 +Each commutation is supported as well. 5.258 + 5.259 +For short distances, the result is very accurate (i.e. respects the dimensions 5.260 +of the WGS-84 spheroid). For longer distances in the order of magnitude of 5.261 +earth's radius or greater, the value is only approximate (but the error is 5.262 +still less than 0.2% as long as no polygons with very long edges are involved). 5.263 + 5.264 +The functions `distance(epoint, epoint)` and `distance(ecluster, epoint)` can 5.265 +be used as an alias for this operator. 5.266 + 5.267 +Note: In case of radial searches with a fixed radius, this operator should 5.268 +not be used. Instead, an `ecircle` should be created and used in combination 5.269 +with the overlap operator (`&&`). Alternatively, the functions 5.270 +`distance_within(epoint, epoint, float8)` or `distance_within(ecluster, epoint, 5.271 +float8)` can be used for fixed-radius searches. 5.272 + 5.273 +### 4. Functions 5.274 + 5.275 +#### `center(circle)` 5.276 + 5.277 +Returns the center of an `ecircle` as an `epoint`. 5.278 + 5.279 +#### `distance(epoint, epoint)` 5.280 + 5.281 +Calculates the distance between two `epoint` datums in meters. This function is 5.282 +an alias for the distance operator `<->`. 5.283 + 5.284 +Note: In case of radial searches with a fixed radius, this function should not be 5.285 +used. Use `distance_within(epoint, epoint, float8)` instead. 5.286 + 5.287 +#### `distance(ecluster, epoint)` 5.288 + 5.289 +Calculates the distance from an `ecluster` to an `epoint` in meters. This 5.290 +function is an alias for the distance operator `<->`. 5.291 + 5.292 +Note: In case of radial searches with a fixed radius, this function should not be 5.293 +used. Use `distance_within(epoint, epoint, float8)` instead. 5.294 + 5.295 +#### `distance_within(`variable `epoint,` fixed `epoint,` radius `float8)` 5.296 + 5.297 +Checks if the distance between two `epoint` datums is not greater than a given 5.298 +value (search radius). 5.299 + 5.300 +Note: In case of radial searches with a fixed radius, the first argument must 5.301 +be used for the table column, while the second argument must be used for the 5.302 +search center. Otherwise an existing index cannot be used. 5.303 + 5.304 +#### `distance_within(`variable `ecluster,` fixed `epoint,` radius `float8)` 5.305 + 5.306 +Checks if the distance from an `ecluster` to an `epoint` is not greater than a 5.307 +given value (search radius). 5.308 + 5.309 +#### `ebox(`latmin `float8,` latmax `float8,` lonmin `float8,` lonmax `float8)` 5.310 + 5.311 +Creates a new `ebox` with the given boundaries. 5.312 +See "1. Types", subsection `ebox` for details. 5.313 + 5.314 +#### `ebox(epoint, epoint)` 5.315 + 5.316 +Creates a new `ebox`. This function may only be used if the longitude 5.317 +difference is less than or equal to 120 degrees. 5.318 +See "1. Types", subsection `ebox` for details. 5.319 + 5.320 +#### `ecircle(epoint, float8)` 5.321 + 5.322 +Creates an `ecircle` with the given center point and radius. 5.323 + 5.324 +#### `ecircle(`latitude `float8,` longitude `float8,` radius `float8)` 5.325 + 5.326 +Creates an `ecircle` with the given center point and radius. 5.327 + 5.328 +#### `ecluster_concat(ecluster, ecluster)` 5.329 + 5.330 +Combines two clusters to form a new `ecluster` by uniting all entries of both 5.331 +clusters. Note that two overlapping areas of polygons annihilate each other 5.332 +(which may be used to create polygons with holes). 5.333 + 5.334 +#### `ecluster_concat(ecluster[])` 5.335 + 5.336 +Creates a new `ecluster` that unites all entries of all clusters in the passed 5.337 +array. Note that two overlapping areas of polygons annihilate each other (which 5.338 +may be used to create polygons with holes). 5.339 + 5.340 +#### `ecluster_create_multipoint(epoint[])` 5.341 + 5.342 +Creates a new `ecluster` which contains multiple points. 5.343 + 5.344 +#### `ecluster_create_outline(epoint[])` 5.345 + 5.346 +Creates a new `ecluster` that is an outline given by the passed points. 5.347 + 5.348 +#### `ecluster_create_path(epoint[])` 5.349 + 5.350 +Creates a new `ecluster` that is a path given by the passed points. 5.351 + 5.352 +#### `ecluster_create_polygon(epoint[])` 5.353 + 5.354 +Creates a new `ecluster` that is a polygon given by the passed points. 5.355 + 5.356 +#### `ecluster_extract_outlines(ecluster)` 5.357 + 5.358 +Set-returning function that returns the outlines of an `ecluster` as `epoint[]` 5.359 +rows. 5.360 + 5.361 +#### `ecluster_extract_paths(ecluster)` 5.362 + 5.363 +Set-returning function that returns the paths of an `ecluster` as `epoint[]` 5.364 +rows. 5.365 + 5.366 +#### `ecluster_extract_points(ecluster)` 5.367 + 5.368 +Set-returning function that returns the points of an `ecluster` as `epoint` 5.369 +rows. 5.370 + 5.371 +#### `ecluster_extract_polygons(ecluster)` 5.372 + 5.373 +Set-returning function that returns the polygons of an `ecluster` as `epoint[]` 5.374 +rows. 5.375 + 5.376 +#### `empty_ebox`() 5.377 + 5.378 +Returns the empty `ebox`. 5.379 +See "1. Types", subsection `ebox` for details. 5.380 + 5.381 +#### `epoint(`latitude `float8,` longitude `float8)` 5.382 + 5.383 +Returns an `epoint` with the given latitude and longitude. 5.384 + 5.385 +#### `epoint_latlon(`latitude `float8,` longitude `float8)` 5.386 + 5.387 +Alias for `epoint(float8, float8)`. 5.388 + 5.389 +#### `epoint_lonlat(`longitude `float8,` latitude `float8)` 5.390 + 5.391 +Same as `epoint(float8, float8)` but with arguments reversed. 5.392 + 5.393 +#### `GeoJSON_to_epoint(jsonb, text)` 5.394 + 5.395 +Maps a GeoJSON object of type "Point" or "Feature" (which contains a 5.396 +"Point") to an `epoint` datum. For any other JSON objects, NULL is returned. 5.397 + 5.398 +The second parameter (which defaults to `epoint_lonlat`) may be set to a name 5.399 +of a conversion function that transforms two coordinates (two `float8` 5.400 +parameters) to an `epoint`. 5.401 + 5.402 +#### `GeoJSON_to_ecluster(jsonb, text)` 5.403 + 5.404 +Maps a (valid) GeoJSON object to an `ecluster`. Note that this function 5.405 +does not check whether the JSONB object is a valid GeoJSON object. 5.406 + 5.407 +The second parameter (which defaults to `epoint_lonlat`) may be set to a name 5.408 +of a conversion function that transforms two coordinates (two `float8` 5.409 +parameters) to an `epoint`. 5.410 + 5.411 +#### `max_latitude(ebox)` 5.412 + 5.413 +Returns the northern boundary of a given `ebox` in degrees between -90 and +90. 5.414 + 5.415 +#### `max_longitude(ebox)` 5.416 + 5.417 +Returns the eastern boundary of a given `ebox` in degrees between -180 and +180 5.418 +(both inclusive). 5.419 + 5.420 +#### `min_latitude(ebox)` 5.421 + 5.422 +Returns the southern boundary of a given `ebox` in degrees between -90 and +90. 5.423 + 5.424 +#### `min_longitude(ebox)` 5.425 + 5.426 +Returns the western boundary of a given `ebox` in degrees between -180 and +180 5.427 +(both inclusive). 5.428 + 5.429 +#### `latitude(epoint)` 5.430 + 5.431 +Returns the latitude value of an `epoint` in degrees between -90 and +90. 5.432 + 5.433 +#### `longitude(epoint)` 5.434 + 5.435 +Returns the longitude value of an `epoint` in degrees between -180 and +180 5.436 +(both inclusive). 5.437 + 5.438 +#### `radius(ecircle)` 5.439 + 5.440 +Returns the radius of an `ecircle` in meters. 5.441 +
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 6.2 +++ b/create_test_db.data.sql Sun Aug 21 17:43:48 2016 +0200 6.3 @@ -0,0 +1,6 @@ 6.4 + 6.5 +INSERT INTO "test" ("location", "area") SELECT 6.6 + epoint((asin(2*random()-1) / pi()) * 180, (2*random()-1) * 180), 6.7 + ecircle((asin(2*random()-1) / pi()) * 180, (2*random()-1) * 180, -ln(1-random()) * 1000) 6.8 + FROM generate_series(1, 10000); 6.9 +
7.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 7.2 +++ b/create_test_db.schema.sql Sun Aug 21 17:43:48 2016 +0200 7.3 @@ -0,0 +1,11 @@ 7.4 + 7.5 +CREATE EXTENSION latlon; 7.6 + 7.7 +CREATE TABLE "test" ( 7.8 + "id" SERIAL4 PRIMARY KEY, 7.9 + "location" EPOINT NOT NULL, 7.10 + "area" ECIRCLE NOT NULL ); 7.11 + 7.12 +CREATE INDEX "test_location_key" ON "test" USING gist ("location"); 7.13 +CREATE INDEX "test_area_key" ON "test" USING gist ("area"); 7.14 +
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 8.2 +++ b/create_test_db.sh Sun Aug 21 17:43:48 2016 +0200 8.3 @@ -0,0 +1,8 @@ 8.4 +#!/bin/sh 8.5 +dropdb latlon_test --if-exists 8.6 +createdb latlon_test || exit 1 8.7 +psql -v ON_ERROR_STOP=1 -f create_test_db.schema.sql latlon_test || exit 1 8.8 +for i in 1 2 3 4 5 6 7 8 9 10 8.9 +do 8.10 +psql -v ON_ERROR_STOP=1 -f create_test_db.data.sql latlon_test || exit 1 8.11 +done
9.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 9.2 +++ b/latlon--0.1.sql Sun Aug 21 17:43:48 2016 +0200 9.3 @@ -0,0 +1,1205 @@ 9.4 + 9.5 +---------------------------------------- 9.6 +-- forward declarations (shell types) -- 9.7 +---------------------------------------- 9.8 + 9.9 +CREATE TYPE epoint; 9.10 +CREATE TYPE ebox; 9.11 +CREATE TYPE ecircle; 9.12 +CREATE TYPE ecluster; 9.13 + 9.14 + 9.15 +------------------------------------------------------------ 9.16 +-- dummy input/output functions for dummy index key types -- 9.17 +------------------------------------------------------------ 9.18 + 9.19 +CREATE FUNCTION ekey_point_in_dummy(cstring) 9.20 + RETURNS ekey_point 9.21 + LANGUAGE C IMMUTABLE STRICT 9.22 + AS '$libdir/latlon-v0001', 'pgl_notimpl'; 9.23 + 9.24 +CREATE FUNCTION ekey_point_out_dummy(ekey_point) 9.25 + RETURNS cstring 9.26 + LANGUAGE C IMMUTABLE STRICT 9.27 + AS '$libdir/latlon-v0001', 'pgl_notimpl'; 9.28 + 9.29 +CREATE FUNCTION ekey_area_in_dummy(cstring) 9.30 + RETURNS ekey_area 9.31 + LANGUAGE C IMMUTABLE STRICT 9.32 + AS '$libdir/latlon-v0001', 'pgl_notimpl'; 9.33 + 9.34 +CREATE FUNCTION ekey_area_out_dummy(ekey_area) 9.35 + RETURNS cstring 9.36 + LANGUAGE C IMMUTABLE STRICT 9.37 + AS '$libdir/latlon-v0001', 'pgl_notimpl'; 9.38 + 9.39 + 9.40 +-------------------------- 9.41 +-- text input functions -- 9.42 +-------------------------- 9.43 + 9.44 +CREATE FUNCTION epoint_in(cstring) 9.45 + RETURNS epoint 9.46 + LANGUAGE C IMMUTABLE STRICT 9.47 + AS '$libdir/latlon-v0001', 'pgl_epoint_in'; 9.48 + 9.49 +CREATE FUNCTION ebox_in(cstring) 9.50 + RETURNS ebox 9.51 + LANGUAGE C IMMUTABLE STRICT 9.52 + AS '$libdir/latlon-v0001', 'pgl_ebox_in'; 9.53 + 9.54 +CREATE FUNCTION ecircle_in(cstring) 9.55 + RETURNS ecircle 9.56 + LANGUAGE C IMMUTABLE STRICT 9.57 + AS '$libdir/latlon-v0001', 'pgl_ecircle_in'; 9.58 + 9.59 +CREATE FUNCTION ecluster_in(cstring) 9.60 + RETURNS ecluster 9.61 + LANGUAGE C IMMUTABLE STRICT 9.62 + AS '$libdir/latlon-v0001', 'pgl_ecluster_in'; 9.63 + 9.64 + 9.65 +--------------------------- 9.66 +-- text output functions -- 9.67 +--------------------------- 9.68 + 9.69 +CREATE FUNCTION epoint_out(epoint) 9.70 + RETURNS cstring 9.71 + LANGUAGE C IMMUTABLE STRICT 9.72 + AS '$libdir/latlon-v0001', 'pgl_epoint_out'; 9.73 + 9.74 +CREATE FUNCTION ebox_out(ebox) 9.75 + RETURNS cstring 9.76 + LANGUAGE C IMMUTABLE STRICT 9.77 + AS '$libdir/latlon-v0001', 'pgl_ebox_out'; 9.78 + 9.79 +CREATE FUNCTION ecircle_out(ecircle) 9.80 + RETURNS cstring 9.81 + LANGUAGE C IMMUTABLE STRICT 9.82 + AS '$libdir/latlon-v0001', 'pgl_ecircle_out'; 9.83 + 9.84 +CREATE FUNCTION ecluster_out(ecluster) 9.85 + RETURNS cstring 9.86 + LANGUAGE C IMMUTABLE STRICT 9.87 + AS '$libdir/latlon-v0001', 'pgl_ecluster_out'; 9.88 + 9.89 + 9.90 +-------------------------- 9.91 +-- binary I/O functions -- 9.92 +-------------------------- 9.93 + 9.94 +CREATE FUNCTION epoint_recv(internal) 9.95 + RETURNS epoint 9.96 + LANGUAGE C IMMUTABLE STRICT 9.97 + AS '$libdir/latlon-v0001', 'pgl_epoint_recv'; 9.98 + 9.99 +CREATE FUNCTION ebox_recv(internal) 9.100 + RETURNS ebox 9.101 + LANGUAGE C IMMUTABLE STRICT 9.102 + AS '$libdir/latlon-v0001', 'pgl_ebox_recv'; 9.103 + 9.104 +CREATE FUNCTION ecircle_recv(internal) 9.105 + RETURNS ecircle 9.106 + LANGUAGE C IMMUTABLE STRICT 9.107 + AS '$libdir/latlon-v0001', 'pgl_ecircle_recv'; 9.108 + 9.109 +CREATE FUNCTION epoint_send(epoint) 9.110 + RETURNS bytea 9.111 + LANGUAGE C IMMUTABLE STRICT 9.112 + AS '$libdir/latlon-v0001', 'pgl_epoint_send'; 9.113 + 9.114 +CREATE FUNCTION ebox_send(ebox) 9.115 + RETURNS bytea 9.116 + LANGUAGE C IMMUTABLE STRICT 9.117 + AS '$libdir/latlon-v0001', 'pgl_ebox_send'; 9.118 + 9.119 +CREATE FUNCTION ecircle_send(ecircle) 9.120 + RETURNS bytea 9.121 + LANGUAGE C IMMUTABLE STRICT 9.122 + AS '$libdir/latlon-v0001', 'pgl_ecircle_send'; 9.123 + 9.124 + 9.125 +----------------------------------------------- 9.126 +-- type definitions of dummy index key types -- 9.127 +----------------------------------------------- 9.128 + 9.129 +CREATE TYPE ekey_point ( 9.130 + internallength = 8, 9.131 + input = ekey_point_in_dummy, 9.132 + output = ekey_point_out_dummy, 9.133 + alignment = char ); 9.134 + 9.135 +CREATE TYPE ekey_area ( 9.136 + internallength = 9, 9.137 + input = ekey_area_in_dummy, 9.138 + output = ekey_area_out_dummy, 9.139 + alignment = char ); 9.140 + 9.141 + 9.142 +------------------------------------------ 9.143 +-- definitions of geographic data types -- 9.144 +------------------------------------------ 9.145 + 9.146 +CREATE TYPE epoint ( 9.147 + internallength = 16, 9.148 + input = epoint_in, 9.149 + output = epoint_out, 9.150 + receive = epoint_recv, 9.151 + send = epoint_send, 9.152 + alignment = double ); 9.153 + 9.154 +CREATE TYPE ebox ( 9.155 + internallength = 32, 9.156 + input = ebox_in, 9.157 + output = ebox_out, 9.158 + receive = ebox_recv, 9.159 + send = ebox_send, 9.160 + alignment = double ); 9.161 + 9.162 +CREATE TYPE ecircle ( 9.163 + internallength = 24, 9.164 + input = ecircle_in, 9.165 + output = ecircle_out, 9.166 + receive = ecircle_recv, 9.167 + send = ecircle_send, 9.168 + alignment = double ); 9.169 + 9.170 +CREATE TYPE ecluster ( 9.171 + internallength = VARIABLE, 9.172 + input = ecluster_in, 9.173 + output = ecluster_out, 9.174 + alignment = double, 9.175 + storage = external ); 9.176 + 9.177 + 9.178 +-------------------- 9.179 +-- B-tree support -- 9.180 +-------------------- 9.181 + 9.182 +-- begin of B-tree support for epoint 9.183 + 9.184 +CREATE FUNCTION epoint_btree_lt(epoint, epoint) 9.185 + RETURNS boolean 9.186 + LANGUAGE C IMMUTABLE STRICT 9.187 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_lt'; 9.188 + 9.189 +CREATE FUNCTION epoint_btree_le(epoint, epoint) 9.190 + RETURNS boolean 9.191 + LANGUAGE C IMMUTABLE STRICT 9.192 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_le'; 9.193 + 9.194 +CREATE FUNCTION epoint_btree_eq(epoint, epoint) 9.195 + RETURNS boolean 9.196 + LANGUAGE C IMMUTABLE STRICT 9.197 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_eq'; 9.198 + 9.199 +CREATE FUNCTION epoint_btree_ne(epoint, epoint) 9.200 + RETURNS boolean 9.201 + LANGUAGE C IMMUTABLE STRICT 9.202 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_ne'; 9.203 + 9.204 +CREATE FUNCTION epoint_btree_ge(epoint, epoint) 9.205 + RETURNS boolean 9.206 + LANGUAGE C IMMUTABLE STRICT 9.207 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_ge'; 9.208 + 9.209 +CREATE FUNCTION epoint_btree_gt(epoint, epoint) 9.210 + RETURNS boolean 9.211 + LANGUAGE C IMMUTABLE STRICT 9.212 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_gt'; 9.213 + 9.214 +CREATE OPERATOR <<< ( 9.215 + leftarg = epoint, 9.216 + rightarg = epoint, 9.217 + procedure = epoint_btree_lt, 9.218 + commutator = >>>, 9.219 + negator = >>>=, 9.220 + restrict = scalarltsel, 9.221 + join = scalarltjoinsel 9.222 +); 9.223 + 9.224 +CREATE OPERATOR <<<= ( 9.225 + leftarg = epoint, 9.226 + rightarg = epoint, 9.227 + procedure = epoint_btree_le, 9.228 + commutator = >>>=, 9.229 + negator = >>>, 9.230 + restrict = scalarltsel, 9.231 + join = scalarltjoinsel 9.232 +); 9.233 + 9.234 +CREATE OPERATOR = ( 9.235 + leftarg = epoint, 9.236 + rightarg = epoint, 9.237 + procedure = epoint_btree_eq, 9.238 + commutator = =, 9.239 + negator = <>, 9.240 + restrict = eqsel, 9.241 + join = eqjoinsel, 9.242 + merges 9.243 +); 9.244 + 9.245 +CREATE OPERATOR <> ( 9.246 + leftarg = epoint, 9.247 + rightarg = epoint, 9.248 + procedure = epoint_btree_eq, 9.249 + commutator = <>, 9.250 + negator = =, 9.251 + restrict = neqsel, 9.252 + join = neqjoinsel 9.253 +); 9.254 + 9.255 +CREATE OPERATOR >>>= ( 9.256 + leftarg = epoint, 9.257 + rightarg = epoint, 9.258 + procedure = epoint_btree_ge, 9.259 + commutator = <<<=, 9.260 + negator = <<<, 9.261 + restrict = scalargtsel, 9.262 + join = scalargtjoinsel 9.263 +); 9.264 + 9.265 +CREATE OPERATOR >>> ( 9.266 + leftarg = epoint, 9.267 + rightarg = epoint, 9.268 + procedure = epoint_btree_gt, 9.269 + commutator = <<<, 9.270 + negator = <<<=, 9.271 + restrict = scalargtsel, 9.272 + join = scalargtjoinsel 9.273 +); 9.274 + 9.275 +CREATE FUNCTION epoint_btree_cmp(epoint, epoint) 9.276 + RETURNS int4 9.277 + LANGUAGE C IMMUTABLE STRICT 9.278 + AS '$libdir/latlon-v0001', 'pgl_btree_epoint_cmp'; 9.279 + 9.280 +CREATE OPERATOR CLASS epoint_btree_ops 9.281 + DEFAULT FOR TYPE epoint USING btree AS 9.282 + OPERATOR 1 <<< , 9.283 + OPERATOR 2 <<<= , 9.284 + OPERATOR 3 = , 9.285 + OPERATOR 4 >>>= , 9.286 + OPERATOR 5 >>> , 9.287 + FUNCTION 1 epoint_btree_cmp(epoint, epoint); 9.288 + 9.289 +-- end of B-tree support for epoint 9.290 + 9.291 +-- begin of B-tree support for ebox 9.292 + 9.293 +CREATE FUNCTION ebox_btree_lt(ebox, ebox) 9.294 + RETURNS boolean 9.295 + LANGUAGE C IMMUTABLE STRICT 9.296 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_lt'; 9.297 + 9.298 +CREATE FUNCTION ebox_btree_le(ebox, ebox) 9.299 + RETURNS boolean 9.300 + LANGUAGE C IMMUTABLE STRICT 9.301 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_le'; 9.302 + 9.303 +CREATE FUNCTION ebox_btree_eq(ebox, ebox) 9.304 + RETURNS boolean 9.305 + LANGUAGE C IMMUTABLE STRICT 9.306 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_eq'; 9.307 + 9.308 +CREATE FUNCTION ebox_btree_ne(ebox, ebox) 9.309 + RETURNS boolean 9.310 + LANGUAGE C IMMUTABLE STRICT 9.311 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_ne'; 9.312 + 9.313 +CREATE FUNCTION ebox_btree_ge(ebox, ebox) 9.314 + RETURNS boolean 9.315 + LANGUAGE C IMMUTABLE STRICT 9.316 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_ge'; 9.317 + 9.318 +CREATE FUNCTION ebox_btree_gt(ebox, ebox) 9.319 + RETURNS boolean 9.320 + LANGUAGE C IMMUTABLE STRICT 9.321 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_gt'; 9.322 + 9.323 +CREATE OPERATOR <<< ( 9.324 + leftarg = ebox, 9.325 + rightarg = ebox, 9.326 + procedure = ebox_btree_lt, 9.327 + commutator = >>>, 9.328 + negator = >>>=, 9.329 + restrict = scalarltsel, 9.330 + join = scalarltjoinsel 9.331 +); 9.332 + 9.333 +CREATE OPERATOR <<<= ( 9.334 + leftarg = ebox, 9.335 + rightarg = ebox, 9.336 + procedure = ebox_btree_le, 9.337 + commutator = >>>=, 9.338 + negator = >>>, 9.339 + restrict = scalarltsel, 9.340 + join = scalarltjoinsel 9.341 +); 9.342 + 9.343 +CREATE OPERATOR = ( 9.344 + leftarg = ebox, 9.345 + rightarg = ebox, 9.346 + procedure = ebox_btree_eq, 9.347 + commutator = =, 9.348 + negator = <>, 9.349 + restrict = eqsel, 9.350 + join = eqjoinsel, 9.351 + merges 9.352 +); 9.353 + 9.354 +CREATE OPERATOR <> ( 9.355 + leftarg = ebox, 9.356 + rightarg = ebox, 9.357 + procedure = ebox_btree_eq, 9.358 + commutator = <>, 9.359 + negator = =, 9.360 + restrict = neqsel, 9.361 + join = neqjoinsel 9.362 +); 9.363 + 9.364 +CREATE OPERATOR >>>= ( 9.365 + leftarg = ebox, 9.366 + rightarg = ebox, 9.367 + procedure = ebox_btree_ge, 9.368 + commutator = <<<=, 9.369 + negator = <<<, 9.370 + restrict = scalargtsel, 9.371 + join = scalargtjoinsel 9.372 +); 9.373 + 9.374 +CREATE OPERATOR >>> ( 9.375 + leftarg = ebox, 9.376 + rightarg = ebox, 9.377 + procedure = ebox_btree_gt, 9.378 + commutator = <<<, 9.379 + negator = <<<=, 9.380 + restrict = scalargtsel, 9.381 + join = scalargtjoinsel 9.382 +); 9.383 + 9.384 +CREATE FUNCTION ebox_btree_cmp(ebox, ebox) 9.385 + RETURNS int4 9.386 + LANGUAGE C IMMUTABLE STRICT 9.387 + AS '$libdir/latlon-v0001', 'pgl_btree_ebox_cmp'; 9.388 + 9.389 +CREATE OPERATOR CLASS ebox_btree_ops 9.390 + DEFAULT FOR TYPE ebox USING btree AS 9.391 + OPERATOR 1 <<< , 9.392 + OPERATOR 2 <<<= , 9.393 + OPERATOR 3 = , 9.394 + OPERATOR 4 >>>= , 9.395 + OPERATOR 5 >>> , 9.396 + FUNCTION 1 ebox_btree_cmp(ebox, ebox); 9.397 + 9.398 +-- end of B-tree support for ebox 9.399 + 9.400 +-- begin of B-tree support for ecircle 9.401 + 9.402 +CREATE FUNCTION ecircle_btree_lt(ecircle, ecircle) 9.403 + RETURNS boolean 9.404 + LANGUAGE C IMMUTABLE STRICT 9.405 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_lt'; 9.406 + 9.407 +CREATE FUNCTION ecircle_btree_le(ecircle, ecircle) 9.408 + RETURNS boolean 9.409 + LANGUAGE C IMMUTABLE STRICT 9.410 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_le'; 9.411 + 9.412 +CREATE FUNCTION ecircle_btree_eq(ecircle, ecircle) 9.413 + RETURNS boolean 9.414 + LANGUAGE C IMMUTABLE STRICT 9.415 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_eq'; 9.416 + 9.417 +CREATE FUNCTION ecircle_btree_ne(ecircle, ecircle) 9.418 + RETURNS boolean 9.419 + LANGUAGE C IMMUTABLE STRICT 9.420 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_ne'; 9.421 + 9.422 +CREATE FUNCTION ecircle_btree_ge(ecircle, ecircle) 9.423 + RETURNS boolean 9.424 + LANGUAGE C IMMUTABLE STRICT 9.425 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_ge'; 9.426 + 9.427 +CREATE FUNCTION ecircle_btree_gt(ecircle, ecircle) 9.428 + RETURNS boolean 9.429 + LANGUAGE C IMMUTABLE STRICT 9.430 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_gt'; 9.431 + 9.432 +CREATE OPERATOR <<< ( 9.433 + leftarg = ecircle, 9.434 + rightarg = ecircle, 9.435 + procedure = ecircle_btree_lt, 9.436 + commutator = >>>, 9.437 + negator = >>>=, 9.438 + restrict = scalarltsel, 9.439 + join = scalarltjoinsel 9.440 +); 9.441 + 9.442 +CREATE OPERATOR <<<= ( 9.443 + leftarg = ecircle, 9.444 + rightarg = ecircle, 9.445 + procedure = ecircle_btree_le, 9.446 + commutator = >>>=, 9.447 + negator = >>>, 9.448 + restrict = scalarltsel, 9.449 + join = scalarltjoinsel 9.450 +); 9.451 + 9.452 +CREATE OPERATOR = ( 9.453 + leftarg = ecircle, 9.454 + rightarg = ecircle, 9.455 + procedure = ecircle_btree_eq, 9.456 + commutator = =, 9.457 + negator = <>, 9.458 + restrict = eqsel, 9.459 + join = eqjoinsel, 9.460 + merges 9.461 +); 9.462 + 9.463 +CREATE OPERATOR <> ( 9.464 + leftarg = ecircle, 9.465 + rightarg = ecircle, 9.466 + procedure = ecircle_btree_eq, 9.467 + commutator = <>, 9.468 + negator = =, 9.469 + restrict = neqsel, 9.470 + join = neqjoinsel 9.471 +); 9.472 + 9.473 +CREATE OPERATOR >>>= ( 9.474 + leftarg = ecircle, 9.475 + rightarg = ecircle, 9.476 + procedure = ecircle_btree_ge, 9.477 + commutator = <<<=, 9.478 + negator = <<<, 9.479 + restrict = scalargtsel, 9.480 + join = scalargtjoinsel 9.481 +); 9.482 + 9.483 +CREATE OPERATOR >>> ( 9.484 + leftarg = ecircle, 9.485 + rightarg = ecircle, 9.486 + procedure = ecircle_btree_gt, 9.487 + commutator = <<<, 9.488 + negator = <<<=, 9.489 + restrict = scalargtsel, 9.490 + join = scalargtjoinsel 9.491 +); 9.492 + 9.493 +CREATE FUNCTION ecircle_btree_cmp(ecircle, ecircle) 9.494 + RETURNS int4 9.495 + LANGUAGE C IMMUTABLE STRICT 9.496 + AS '$libdir/latlon-v0001', 'pgl_btree_ecircle_cmp'; 9.497 + 9.498 +CREATE OPERATOR CLASS ecircle_btree_ops 9.499 + DEFAULT FOR TYPE ecircle USING btree AS 9.500 + OPERATOR 1 <<< , 9.501 + OPERATOR 2 <<<= , 9.502 + OPERATOR 3 = , 9.503 + OPERATOR 4 >>>= , 9.504 + OPERATOR 5 >>> , 9.505 + FUNCTION 1 ecircle_btree_cmp(ecircle, ecircle); 9.506 + 9.507 +-- end of B-tree support for ecircle 9.508 + 9.509 + 9.510 +---------------- 9.511 +-- type casts -- 9.512 +---------------- 9.513 + 9.514 +CREATE FUNCTION cast_epoint_to_ebox(epoint) 9.515 + RETURNS ebox 9.516 + LANGUAGE C IMMUTABLE STRICT 9.517 + AS '$libdir/latlon-v0001', 'pgl_epoint_to_ebox'; 9.518 + 9.519 +CREATE CAST (epoint AS ebox) WITH FUNCTION cast_epoint_to_ebox(epoint); 9.520 + 9.521 +CREATE FUNCTION cast_epoint_to_ecircle(epoint) 9.522 + RETURNS ecircle 9.523 + LANGUAGE C IMMUTABLE STRICT 9.524 + AS '$libdir/latlon-v0001', 'pgl_epoint_to_ecircle'; 9.525 + 9.526 +CREATE CAST (epoint AS ecircle) WITH FUNCTION cast_epoint_to_ecircle(epoint); 9.527 + 9.528 +CREATE FUNCTION cast_epoint_to_ecluster(epoint) 9.529 + RETURNS ecluster 9.530 + LANGUAGE C IMMUTABLE STRICT 9.531 + AS '$libdir/latlon-v0001', 'pgl_epoint_to_ecluster'; 9.532 + 9.533 +CREATE CAST (epoint AS ecluster) WITH FUNCTION cast_epoint_to_ecluster(epoint); 9.534 + 9.535 +CREATE FUNCTION cast_ebox_to_ecluster(ebox) 9.536 + RETURNS ecluster 9.537 + LANGUAGE C IMMUTABLE STRICT 9.538 + AS '$libdir/latlon-v0001', 'pgl_ebox_to_ecluster'; 9.539 + 9.540 +CREATE CAST (ebox AS ecluster) WITH FUNCTION cast_ebox_to_ecluster(ebox); 9.541 + 9.542 + 9.543 +--------------------------- 9.544 +-- constructor functions -- 9.545 +--------------------------- 9.546 + 9.547 +CREATE FUNCTION epoint(float8, float8) 9.548 + RETURNS epoint 9.549 + LANGUAGE C IMMUTABLE STRICT 9.550 + AS '$libdir/latlon-v0001', 'pgl_create_epoint'; 9.551 + 9.552 +CREATE FUNCTION epoint_latlon(float8, float8) 9.553 + RETURNS epoint 9.554 + LANGUAGE SQL IMMUTABLE STRICT AS $$ 9.555 + SELECT epoint($1, $2) 9.556 + $$; 9.557 + 9.558 +CREATE FUNCTION epoint_lonlat(float8, float8) 9.559 + RETURNS epoint 9.560 + LANGUAGE SQL IMMUTABLE STRICT AS $$ 9.561 + SELECT epoint($2, $1) 9.562 + $$; 9.563 + 9.564 +CREATE FUNCTION empty_ebox() 9.565 + RETURNS ebox 9.566 + LANGUAGE C IMMUTABLE STRICT 9.567 + AS '$libdir/latlon-v0001', 'pgl_create_empty_ebox'; 9.568 + 9.569 +CREATE FUNCTION ebox(float8, float8, float8, float8) 9.570 + RETURNS ebox 9.571 + LANGUAGE C IMMUTABLE STRICT 9.572 + AS '$libdir/latlon-v0001', 'pgl_create_ebox'; 9.573 + 9.574 +CREATE FUNCTION ebox(epoint, epoint) 9.575 + RETURNS ebox 9.576 + LANGUAGE C IMMUTABLE STRICT 9.577 + AS '$libdir/latlon-v0001', 'pgl_create_ebox_from_epoints'; 9.578 + 9.579 +CREATE FUNCTION ecircle(float8, float8, float8) 9.580 + RETURNS ecircle 9.581 + LANGUAGE C IMMUTABLE STRICT 9.582 + AS '$libdir/latlon-v0001', 'pgl_create_ecircle'; 9.583 + 9.584 +CREATE FUNCTION ecircle(epoint, float8) 9.585 + RETURNS ecircle 9.586 + LANGUAGE C IMMUTABLE STRICT 9.587 + AS '$libdir/latlon-v0001', 'pgl_create_ecircle_from_epoint'; 9.588 + 9.589 +CREATE FUNCTION ecluster_concat(ecluster[]) 9.590 + RETURNS ecluster 9.591 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.592 + SELECT array_to_string($1, ' ')::ecluster 9.593 + $$; 9.594 + 9.595 +CREATE FUNCTION ecluster_concat(ecluster, ecluster) 9.596 + RETURNS ecluster 9.597 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.598 + SELECT ($1::text || ' ' || $2::text)::ecluster 9.599 + $$; 9.600 + 9.601 +CREATE FUNCTION ecluster_create_multipoint(epoint[]) 9.602 + RETURNS ecluster 9.603 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.604 + SELECT 9.605 + array_to_string(array_agg('point (' || unnest || ')'), ' ')::ecluster 9.606 + FROM unnest($1) 9.607 + $$; 9.608 + 9.609 +CREATE FUNCTION ecluster_create_path(epoint[]) 9.610 + RETURNS ecluster 9.611 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.612 + SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE 9.613 + ('path (' || array_to_string($1, ' ') || ')')::ecluster 9.614 + END 9.615 + FROM array_to_string($1, ' ') AS "str" 9.616 + $$; 9.617 + 9.618 +CREATE FUNCTION ecluster_create_outline(epoint[]) 9.619 + RETURNS ecluster 9.620 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.621 + SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE 9.622 + ('outline (' || array_to_string($1, ' ') || ')')::ecluster 9.623 + END 9.624 + FROM array_to_string($1, ' ') AS "str" 9.625 + $$; 9.626 + 9.627 +CREATE FUNCTION ecluster_create_polygon(epoint[]) 9.628 + RETURNS ecluster 9.629 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.630 + SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE 9.631 + ('polygon (' || array_to_string($1, ' ') || ')')::ecluster 9.632 + END 9.633 + FROM array_to_string($1, ' ') AS "str" 9.634 + $$; 9.635 + 9.636 + 9.637 +---------------------- 9.638 +-- getter functions -- 9.639 +---------------------- 9.640 + 9.641 +CREATE FUNCTION latitude(epoint) 9.642 + RETURNS float8 9.643 + LANGUAGE C IMMUTABLE STRICT 9.644 + AS '$libdir/latlon-v0001', 'pgl_epoint_lat'; 9.645 + 9.646 +CREATE FUNCTION longitude(epoint) 9.647 + RETURNS float8 9.648 + LANGUAGE C IMMUTABLE STRICT 9.649 + AS '$libdir/latlon-v0001', 'pgl_epoint_lon'; 9.650 + 9.651 +CREATE FUNCTION min_latitude(ebox) 9.652 + RETURNS float8 9.653 + LANGUAGE C IMMUTABLE STRICT 9.654 + AS '$libdir/latlon-v0001', 'pgl_ebox_lat_min'; 9.655 + 9.656 +CREATE FUNCTION max_latitude(ebox) 9.657 + RETURNS float8 9.658 + LANGUAGE C IMMUTABLE STRICT 9.659 + AS '$libdir/latlon-v0001', 'pgl_ebox_lat_max'; 9.660 + 9.661 +CREATE FUNCTION min_longitude(ebox) 9.662 + RETURNS float8 9.663 + LANGUAGE C IMMUTABLE STRICT 9.664 + AS '$libdir/latlon-v0001', 'pgl_ebox_lon_min'; 9.665 + 9.666 +CREATE FUNCTION max_longitude(ebox) 9.667 + RETURNS float8 9.668 + LANGUAGE C IMMUTABLE STRICT 9.669 + AS '$libdir/latlon-v0001', 'pgl_ebox_lon_max'; 9.670 + 9.671 +CREATE FUNCTION center(ecircle) 9.672 + RETURNS epoint 9.673 + LANGUAGE C IMMUTABLE STRICT 9.674 + AS '$libdir/latlon-v0001', 'pgl_ecircle_center'; 9.675 + 9.676 +CREATE FUNCTION radius(ecircle) 9.677 + RETURNS float8 9.678 + LANGUAGE C IMMUTABLE STRICT 9.679 + AS '$libdir/latlon-v0001', 'pgl_ecircle_radius'; 9.680 + 9.681 +CREATE FUNCTION ecluster_extract_points(ecluster) 9.682 + RETURNS SETOF epoint 9.683 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.684 + SELECT "match"[2]::epoint 9.685 + FROM regexp_matches($1::text, e'(^| )point \\(([^)]+)\\)', 'g') AS "match" 9.686 + $$; 9.687 + 9.688 +CREATE FUNCTION ecluster_extract_paths(ecluster) 9.689 + RETURNS SETOF epoint[] 9.690 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.691 + SELECT ( 9.692 + SELECT array_agg("m2"[1]::epoint) 9.693 + FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2" 9.694 + ) 9.695 + FROM regexp_matches($1::text, e'(^| )path \\(([^)]+)\\)', 'g') AS "m1" 9.696 + $$; 9.697 + 9.698 +CREATE FUNCTION ecluster_extract_outlines(ecluster) 9.699 + RETURNS SETOF epoint[] 9.700 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.701 + SELECT ( 9.702 + SELECT array_agg("m2"[1]::epoint) 9.703 + FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2" 9.704 + ) 9.705 + FROM regexp_matches($1::text, e'(^| )outline \\(([^)]+)\\)', 'g') AS "m1" 9.706 + $$; 9.707 + 9.708 +CREATE FUNCTION ecluster_extract_polygons(ecluster) 9.709 + RETURNS SETOF epoint[] 9.710 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.711 + SELECT ( 9.712 + SELECT array_agg("m2"[1]::epoint) 9.713 + FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2" 9.714 + ) 9.715 + FROM regexp_matches($1::text, e'(^| )polygon \\(([^)]+)\\)', 'g') AS "m1" 9.716 + $$; 9.717 + 9.718 + 9.719 +--------------- 9.720 +-- operators -- 9.721 +--------------- 9.722 + 9.723 +CREATE FUNCTION epoint_ebox_overlap_proc(epoint, ebox) 9.724 + RETURNS boolean 9.725 + LANGUAGE C IMMUTABLE STRICT 9.726 + AS '$libdir/latlon-v0001', 'pgl_epoint_ebox_overlap'; 9.727 + 9.728 +CREATE FUNCTION epoint_ecircle_overlap_proc(epoint, ecircle) 9.729 + RETURNS boolean 9.730 + LANGUAGE C IMMUTABLE STRICT 9.731 + AS '$libdir/latlon-v0001', 'pgl_epoint_ecircle_overlap'; 9.732 + 9.733 +CREATE FUNCTION epoint_ecluster_overlap_proc(epoint, ecluster) 9.734 + RETURNS boolean 9.735 + LANGUAGE C IMMUTABLE STRICT 9.736 + AS '$libdir/latlon-v0001', 'pgl_epoint_ecluster_overlap'; 9.737 + 9.738 +CREATE FUNCTION ebox_overlap_proc(ebox, ebox) 9.739 + RETURNS boolean 9.740 + LANGUAGE C IMMUTABLE STRICT 9.741 + AS '$libdir/latlon-v0001', 'pgl_ebox_overlap'; 9.742 + 9.743 +CREATE FUNCTION ecircle_overlap_proc(ecircle, ecircle) 9.744 + RETURNS boolean 9.745 + LANGUAGE C IMMUTABLE STRICT 9.746 + AS '$libdir/latlon-v0001', 'pgl_ecircle_overlap'; 9.747 + 9.748 +CREATE FUNCTION ecircle_ecluster_overlap_proc(ecircle, ecluster) 9.749 + RETURNS boolean 9.750 + LANGUAGE C IMMUTABLE STRICT 9.751 + AS '$libdir/latlon-v0001', 'pgl_ecircle_ecluster_overlap'; 9.752 + 9.753 +CREATE FUNCTION epoint_distance_proc(epoint, epoint) 9.754 + RETURNS float8 9.755 + LANGUAGE C IMMUTABLE STRICT 9.756 + AS '$libdir/latlon-v0001', 'pgl_epoint_distance'; 9.757 + 9.758 +CREATE FUNCTION epoint_ecircle_distance_proc(epoint, ecircle) 9.759 + RETURNS float8 9.760 + LANGUAGE C IMMUTABLE STRICT 9.761 + AS '$libdir/latlon-v0001', 'pgl_epoint_ecircle_distance'; 9.762 + 9.763 +CREATE FUNCTION epoint_ecluster_distance_proc(epoint, ecluster) 9.764 + RETURNS float8 9.765 + LANGUAGE C IMMUTABLE STRICT 9.766 + AS '$libdir/latlon-v0001', 'pgl_epoint_ecluster_distance'; 9.767 + 9.768 +CREATE FUNCTION ecircle_distance_proc(ecircle, ecircle) 9.769 + RETURNS float8 9.770 + LANGUAGE C IMMUTABLE STRICT 9.771 + AS '$libdir/latlon-v0001', 'pgl_ecircle_distance'; 9.772 + 9.773 +CREATE FUNCTION ecircle_ecluster_distance_proc(ecircle, ecluster) 9.774 + RETURNS float8 9.775 + LANGUAGE C IMMUTABLE STRICT 9.776 + AS '$libdir/latlon-v0001', 'pgl_ecircle_ecluster_distance'; 9.777 + 9.778 +CREATE OPERATOR && ( 9.779 + leftarg = epoint, 9.780 + rightarg = ebox, 9.781 + procedure = epoint_ebox_overlap_proc, 9.782 + commutator = &&, 9.783 + restrict = areasel, 9.784 + join = areajoinsel 9.785 +); 9.786 + 9.787 +CREATE FUNCTION epoint_ebox_overlap_commutator(ebox, epoint) 9.788 + RETURNS boolean 9.789 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 9.790 + 9.791 +CREATE OPERATOR && ( 9.792 + leftarg = ebox, 9.793 + rightarg = epoint, 9.794 + procedure = epoint_ebox_overlap_commutator, 9.795 + commutator = &&, 9.796 + restrict = areasel, 9.797 + join = areajoinsel 9.798 +); 9.799 + 9.800 +CREATE OPERATOR && ( 9.801 + leftarg = epoint, 9.802 + rightarg = ecircle, 9.803 + procedure = epoint_ecircle_overlap_proc, 9.804 + commutator = &&, 9.805 + restrict = areasel, 9.806 + join = areajoinsel 9.807 +); 9.808 + 9.809 +CREATE FUNCTION epoint_ecircle_overlap_commutator(ecircle, epoint) 9.810 + RETURNS boolean 9.811 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 9.812 + 9.813 +CREATE OPERATOR && ( 9.814 + leftarg = ecircle, 9.815 + rightarg = epoint, 9.816 + procedure = epoint_ecircle_overlap_commutator, 9.817 + commutator = &&, 9.818 + restrict = areasel, 9.819 + join = areajoinsel 9.820 +); 9.821 + 9.822 +CREATE OPERATOR && ( 9.823 + leftarg = epoint, 9.824 + rightarg = ecluster, 9.825 + procedure = epoint_ecluster_overlap_proc, 9.826 + commutator = &&, 9.827 + restrict = areasel, 9.828 + join = areajoinsel 9.829 +); 9.830 + 9.831 +CREATE FUNCTION epoint_ecluster_overlap_commutator(ecluster, epoint) 9.832 + RETURNS boolean 9.833 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 9.834 + 9.835 +CREATE OPERATOR && ( 9.836 + leftarg = ecluster, 9.837 + rightarg = epoint, 9.838 + procedure = epoint_ecluster_overlap_commutator, 9.839 + commutator = &&, 9.840 + restrict = areasel, 9.841 + join = areajoinsel 9.842 +); 9.843 + 9.844 +CREATE OPERATOR && ( 9.845 + leftarg = ebox, 9.846 + rightarg = ebox, 9.847 + procedure = ebox_overlap_proc, 9.848 + commutator = &&, 9.849 + restrict = areasel, 9.850 + join = areajoinsel 9.851 +); 9.852 + 9.853 +CREATE OPERATOR && ( 9.854 + leftarg = ecircle, 9.855 + rightarg = ecircle, 9.856 + procedure = ecircle_overlap_proc, 9.857 + commutator = &&, 9.858 + restrict = areasel, 9.859 + join = areajoinsel 9.860 +); 9.861 + 9.862 +CREATE OPERATOR && ( 9.863 + leftarg = ecircle, 9.864 + rightarg = ecluster, 9.865 + procedure = ecircle_ecluster_overlap_proc, 9.866 + commutator = &&, 9.867 + restrict = areasel, 9.868 + join = areajoinsel 9.869 +); 9.870 + 9.871 +CREATE FUNCTION ecircle_ecluster_overlap_commutator(ecluster, ecircle) 9.872 + RETURNS boolean 9.873 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1'; 9.874 + 9.875 +CREATE OPERATOR && ( 9.876 + leftarg = ecluster, 9.877 + rightarg = ecircle, 9.878 + procedure = ecircle_ecluster_overlap_commutator, 9.879 + commutator = &&, 9.880 + restrict = areasel, 9.881 + join = areajoinsel 9.882 +); 9.883 + 9.884 +CREATE OPERATOR <-> ( 9.885 + leftarg = epoint, 9.886 + rightarg = epoint, 9.887 + procedure = epoint_distance_proc, 9.888 + commutator = <-> 9.889 +); 9.890 + 9.891 +CREATE OPERATOR <-> ( 9.892 + leftarg = epoint, 9.893 + rightarg = ecircle, 9.894 + procedure = epoint_ecircle_distance_proc, 9.895 + commutator = <-> 9.896 +); 9.897 + 9.898 +CREATE FUNCTION epoint_ecircle_distance_commutator(ecircle, epoint) 9.899 + RETURNS float8 9.900 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1'; 9.901 + 9.902 +CREATE OPERATOR <-> ( 9.903 + leftarg = ecircle, 9.904 + rightarg = epoint, 9.905 + procedure = epoint_ecircle_distance_commutator, 9.906 + commutator = <-> 9.907 +); 9.908 + 9.909 +CREATE OPERATOR <-> ( 9.910 + leftarg = epoint, 9.911 + rightarg = ecluster, 9.912 + procedure = epoint_ecluster_distance_proc, 9.913 + commutator = <-> 9.914 +); 9.915 + 9.916 +CREATE FUNCTION epoint_ecluster_distance_commutator(ecluster, epoint) 9.917 + RETURNS float8 9.918 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1'; 9.919 + 9.920 +CREATE OPERATOR <-> ( 9.921 + leftarg = ecluster, 9.922 + rightarg = epoint, 9.923 + procedure = epoint_ecluster_distance_commutator, 9.924 + commutator = <-> 9.925 +); 9.926 + 9.927 +CREATE OPERATOR <-> ( 9.928 + leftarg = ecircle, 9.929 + rightarg = ecircle, 9.930 + procedure = ecircle_distance_proc, 9.931 + commutator = <-> 9.932 +); 9.933 + 9.934 +CREATE OPERATOR <-> ( 9.935 + leftarg = ecircle, 9.936 + rightarg = ecluster, 9.937 + procedure = ecircle_ecluster_distance_proc, 9.938 + commutator = <-> 9.939 +); 9.940 + 9.941 +CREATE FUNCTION ecircle_ecluster_distance_commutator(ecluster, ecircle) 9.942 + RETURNS float8 9.943 + LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1'; 9.944 + 9.945 +CREATE OPERATOR <-> ( 9.946 + leftarg = ecluster, 9.947 + rightarg = ecircle, 9.948 + procedure = ecircle_ecluster_distance_commutator, 9.949 + commutator = <-> 9.950 +); 9.951 + 9.952 + 9.953 +---------------- 9.954 +-- GiST index -- 9.955 +---------------- 9.956 + 9.957 +CREATE FUNCTION pgl_gist_consistent(internal, internal, smallint, oid, internal) 9.958 + RETURNS boolean 9.959 + LANGUAGE C STRICT 9.960 + AS '$libdir/latlon-v0001', 'pgl_gist_consistent'; 9.961 + 9.962 +CREATE FUNCTION pgl_gist_union(internal, internal) 9.963 + RETURNS internal 9.964 + LANGUAGE C STRICT 9.965 + AS '$libdir/latlon-v0001', 'pgl_gist_union'; 9.966 + 9.967 +CREATE FUNCTION pgl_gist_compress_epoint(internal) 9.968 + RETURNS internal 9.969 + LANGUAGE C STRICT 9.970 + AS '$libdir/latlon-v0001', 'pgl_gist_compress_epoint'; 9.971 + 9.972 +CREATE FUNCTION pgl_gist_compress_ecircle(internal) 9.973 + RETURNS internal 9.974 + LANGUAGE C STRICT 9.975 + AS '$libdir/latlon-v0001', 'pgl_gist_compress_ecircle'; 9.976 + 9.977 +CREATE FUNCTION pgl_gist_compress_ecluster(internal) 9.978 + RETURNS internal 9.979 + LANGUAGE C STRICT 9.980 + AS '$libdir/latlon-v0001', 'pgl_gist_compress_ecluster'; 9.981 + 9.982 +CREATE FUNCTION pgl_gist_decompress(internal) 9.983 + RETURNS internal 9.984 + LANGUAGE C STRICT 9.985 + AS '$libdir/latlon-v0001', 'pgl_gist_decompress'; 9.986 + 9.987 +CREATE FUNCTION pgl_gist_penalty(internal, internal, internal) 9.988 + RETURNS internal 9.989 + LANGUAGE C STRICT 9.990 + AS '$libdir/latlon-v0001', 'pgl_gist_penalty'; 9.991 + 9.992 +CREATE FUNCTION pgl_gist_picksplit(internal, internal) 9.993 + RETURNS internal 9.994 + LANGUAGE C STRICT 9.995 + AS '$libdir/latlon-v0001', 'pgl_gist_picksplit'; 9.996 + 9.997 +CREATE FUNCTION pgl_gist_same(internal, internal, internal) 9.998 + RETURNS internal 9.999 + LANGUAGE C STRICT 9.1000 + AS '$libdir/latlon-v0001', 'pgl_gist_same'; 9.1001 + 9.1002 +CREATE FUNCTION pgl_gist_distance(internal, internal, smallint, oid) 9.1003 + RETURNS internal 9.1004 + LANGUAGE C STRICT 9.1005 + AS '$libdir/latlon-v0001', 'pgl_gist_distance'; 9.1006 + 9.1007 +CREATE OPERATOR CLASS epoint_ops 9.1008 + DEFAULT FOR TYPE epoint USING gist AS 9.1009 + OPERATOR 11 = , 9.1010 + OPERATOR 22 && (epoint, ebox), 9.1011 + OPERATOR 23 && (epoint, ecircle), 9.1012 + OPERATOR 24 && (epoint, ecluster), 9.1013 + OPERATOR 31 <-> (epoint, epoint) FOR ORDER BY float_ops, 9.1014 + OPERATOR 33 <-> (epoint, ecircle) FOR ORDER BY float_ops, 9.1015 + OPERATOR 34 <-> (epoint, ecluster) FOR ORDER BY float_ops, 9.1016 + FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal), 9.1017 + FUNCTION 2 pgl_gist_union(internal, internal), 9.1018 + FUNCTION 3 pgl_gist_compress_epoint(internal), 9.1019 + FUNCTION 4 pgl_gist_decompress(internal), 9.1020 + FUNCTION 5 pgl_gist_penalty(internal, internal, internal), 9.1021 + FUNCTION 6 pgl_gist_picksplit(internal, internal), 9.1022 + FUNCTION 7 pgl_gist_same(internal, internal, internal), 9.1023 + FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid), 9.1024 + STORAGE ekey_point; 9.1025 + 9.1026 +CREATE OPERATOR CLASS ecircle_ops 9.1027 + DEFAULT FOR TYPE ecircle USING gist AS 9.1028 + OPERATOR 13 = , 9.1029 + OPERATOR 21 && (ecircle, epoint), 9.1030 + OPERATOR 23 && (ecircle, ecircle), 9.1031 + OPERATOR 24 && (ecircle, ecluster), 9.1032 + OPERATOR 31 <-> (ecircle, epoint) FOR ORDER BY float_ops, 9.1033 + OPERATOR 33 <-> (ecircle, ecircle) FOR ORDER BY float_ops, 9.1034 + OPERATOR 34 <-> (ecircle, ecluster) FOR ORDER BY float_ops, 9.1035 + FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal), 9.1036 + FUNCTION 2 pgl_gist_union(internal, internal), 9.1037 + FUNCTION 3 pgl_gist_compress_ecircle(internal), 9.1038 + FUNCTION 4 pgl_gist_decompress(internal), 9.1039 + FUNCTION 5 pgl_gist_penalty(internal, internal, internal), 9.1040 + FUNCTION 6 pgl_gist_picksplit(internal, internal), 9.1041 + FUNCTION 7 pgl_gist_same(internal, internal, internal), 9.1042 + FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid), 9.1043 + STORAGE ekey_area; 9.1044 + 9.1045 +CREATE OPERATOR CLASS ecluster_ops 9.1046 + DEFAULT FOR TYPE ecluster USING gist AS 9.1047 + OPERATOR 21 && (ecluster, epoint), 9.1048 + FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal), 9.1049 + FUNCTION 2 pgl_gist_union(internal, internal), 9.1050 + FUNCTION 3 pgl_gist_compress_ecluster(internal), 9.1051 + FUNCTION 4 pgl_gist_decompress(internal), 9.1052 + FUNCTION 5 pgl_gist_penalty(internal, internal, internal), 9.1053 + FUNCTION 6 pgl_gist_picksplit(internal, internal), 9.1054 + FUNCTION 7 pgl_gist_same(internal, internal, internal), 9.1055 + FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid), 9.1056 + STORAGE ekey_area; 9.1057 + 9.1058 + 9.1059 +--------------------- 9.1060 +-- alias functions -- 9.1061 +--------------------- 9.1062 + 9.1063 +CREATE FUNCTION distance(epoint, epoint) 9.1064 + RETURNS float8 9.1065 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2'; 9.1066 + 9.1067 +CREATE FUNCTION distance(ecluster, epoint) 9.1068 + RETURNS float8 9.1069 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2'; 9.1070 + 9.1071 +CREATE FUNCTION distance_within(epoint, epoint, float8) 9.1072 + RETURNS boolean 9.1073 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)'; 9.1074 + 9.1075 +CREATE FUNCTION distance_within(ecluster, epoint, float8) 9.1076 + RETURNS boolean 9.1077 + LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)'; 9.1078 + 9.1079 + 9.1080 +-------------------------------- 9.1081 +-- other data storage formats -- 9.1082 +-------------------------------- 9.1083 + 9.1084 +CREATE FUNCTION coords_to_epoint(float8, float8, text = 'epoint_lonlat') 9.1085 + RETURNS epoint 9.1086 + LANGUAGE plpgsql IMMUTABLE STRICT AS $$ 9.1087 + DECLARE 9.1088 + "result" epoint; 9.1089 + BEGIN 9.1090 + IF $3 = 'epoint_lonlat' THEN 9.1091 + -- avoid dynamic command execution for better performance 9.1092 + RETURN epoint($2, $1); 9.1093 + END IF; 9.1094 + IF $3 = 'epoint' OR $3 = 'epoint_latlon' THEN 9.1095 + -- avoid dynamic command execution for better performance 9.1096 + RETURN epoint($1, $2); 9.1097 + END IF; 9.1098 + EXECUTE 'SELECT ' || $3 || '($1, $2)' INTO STRICT "result" USING $1, $2; 9.1099 + RETURN "result"; 9.1100 + END; 9.1101 + $$; 9.1102 + 9.1103 +CREATE FUNCTION GeoJSON_to_epoint(jsonb, text = 'epoint_lonlat') 9.1104 + RETURNS epoint 9.1105 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.1106 + SELECT CASE 9.1107 + WHEN $1->>'type' = 'Point' THEN 9.1108 + coords_to_epoint( 9.1109 + ($1->'coordinates'->>1)::float8, 9.1110 + ($1->'coordinates'->>0)::float8, 9.1111 + $2 9.1112 + ) 9.1113 + WHEN $1->>'type' = 'Feature' THEN 9.1114 + GeoJSON_to_epoint($1->'geometry', $2) 9.1115 + ELSE 9.1116 + NULL 9.1117 + END 9.1118 + $$; 9.1119 + 9.1120 +CREATE FUNCTION GeoJSON_to_ecluster(jsonb, text = 'epoint_lonlat') 9.1121 + RETURNS ecluster 9.1122 + LANGUAGE sql IMMUTABLE STRICT AS $$ 9.1123 + SELECT CASE $1->>'type' 9.1124 + WHEN 'Point' THEN 9.1125 + coords_to_epoint( 9.1126 + ($1->'coordinates'->>1)::float8, 9.1127 + ($1->'coordinates'->>0)::float8, 9.1128 + $2 9.1129 + )::ecluster 9.1130 + WHEN 'MultiPoint' THEN 9.1131 + ( SELECT ecluster_create_multipoint(array_agg( 9.1132 + coords_to_epoint( 9.1133 + ("coord"->>1)::float8, 9.1134 + ("coord"->>0)::float8, 9.1135 + $2 9.1136 + ) 9.1137 + )) 9.1138 + FROM jsonb_array_elements($1->'coordinates') AS "coord" 9.1139 + ) 9.1140 + WHEN 'LineString' THEN 9.1141 + ( SELECT ecluster_create_path(array_agg( 9.1142 + coords_to_epoint( 9.1143 + ("coord"->>1)::float8, 9.1144 + ("coord"->>0)::float8, 9.1145 + $2 9.1146 + ) 9.1147 + )) 9.1148 + FROM jsonb_array_elements($1->'coordinates') AS "coord" 9.1149 + ) 9.1150 + WHEN 'MultiLineString' THEN 9.1151 + ( SELECT ecluster_concat(array_agg( 9.1152 + ( SELECT ecluster_create_path(array_agg( 9.1153 + coords_to_epoint( 9.1154 + ("coord"->>1)::float8, 9.1155 + ("coord"->>0)::float8, 9.1156 + $2 9.1157 + ) 9.1158 + )) 9.1159 + FROM jsonb_array_elements("coord_array") AS "coord" 9.1160 + ) 9.1161 + )) 9.1162 + FROM jsonb_array_elements($1->'coordinates') AS "coord_array" 9.1163 + ) 9.1164 + WHEN 'Polygon' THEN 9.1165 + ( SELECT ecluster_concat(array_agg( 9.1166 + ( SELECT ecluster_create_polygon(array_agg( 9.1167 + coords_to_epoint( 9.1168 + ("coord"->>1)::float8, 9.1169 + ("coord"->>0)::float8, 9.1170 + $2 9.1171 + ) 9.1172 + )) 9.1173 + FROM jsonb_array_elements("coord_array") AS "coord" 9.1174 + ) 9.1175 + )) 9.1176 + FROM jsonb_array_elements($1->'coordinates') AS "coord_array" 9.1177 + ) 9.1178 + WHEN 'MultiPolygon' THEN 9.1179 + ( SELECT ecluster_concat(array_agg( 9.1180 + ( SELECT ecluster_concat(array_agg( 9.1181 + ( SELECT ecluster_create_polygon(array_agg( 9.1182 + coords_to_epoint( 9.1183 + ("coord"->>1)::float8, 9.1184 + ("coord"->>0)::float8, 9.1185 + $2 9.1186 + ) 9.1187 + )) 9.1188 + FROM jsonb_array_elements("coord_array") AS "coord" 9.1189 + ) 9.1190 + )) 9.1191 + FROM jsonb_array_elements("coord_array_array") AS "coord_array" 9.1192 + ) 9.1193 + )) 9.1194 + FROM jsonb_array_elements($1->'coordinates') AS "coord_array_array" 9.1195 + ) 9.1196 + WHEN 'Feature' THEN 9.1197 + GeoJSON_to_ecluster($1->'geometry', $2) 9.1198 + WHEN 'FeatureCollection' THEN 9.1199 + ( SELECT ecluster_concat(array_agg( 9.1200 + GeoJSON_to_ecluster("feature", $2) 9.1201 + )) 9.1202 + FROM jsonb_array_elements($1->'features') AS "feature" 9.1203 + ) 9.1204 + ELSE 9.1205 + NULL 9.1206 + END 9.1207 + $$; 9.1208 +
10.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 10.2 +++ b/latlon-v0001.c Sun Aug 21 17:43:48 2016 +0200 10.3 @@ -0,0 +1,2710 @@ 10.4 + 10.5 +/*-------------* 10.6 + * C prelude * 10.7 + *-------------*/ 10.8 + 10.9 +#include "postgres.h" 10.10 +#include "fmgr.h" 10.11 +#include "libpq/pqformat.h" 10.12 +#include "access/gist.h" 10.13 +#include "access/stratnum.h" 10.14 +#include "utils/array.h" 10.15 +#include <math.h> 10.16 + 10.17 +#ifdef PG_MODULE_MAGIC 10.18 +PG_MODULE_MAGIC; 10.19 +#endif 10.20 + 10.21 +#if INT_MAX < 2147483647 10.22 +#error Expected int type to be at least 32 bit wide 10.23 +#endif 10.24 + 10.25 + 10.26 +/*---------------------------------* 10.27 + * distance calculation on earth * 10.28 + * (using WGS-84 spheroid) * 10.29 + *---------------------------------*/ 10.30 + 10.31 +/* WGS-84 spheroid with following parameters: 10.32 + semi-major axis a = 6378137 10.33 + semi-minor axis b = a * (1 - 1/298.257223563) 10.34 + estimated diameter = 2 * (2*a+b)/3 10.35 +*/ 10.36 +#define PGL_SPHEROID_A 6378137.0 /* semi major axis */ 10.37 +#define PGL_SPHEROID_F (1.0/298.257223563) /* flattening */ 10.38 +#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F)) 10.39 +#define PGL_EPS2 ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \ 10.40 + PGL_SPHEROID_B * PGL_SPHEROID_B ) / \ 10.41 + ( PGL_SPHEROID_A * PGL_SPHEROID_A ) ) 10.42 +#define PGL_SUBEPS2 (1.0-PGL_EPS2) 10.43 +#define PGL_DIAMETER ((4.0*PGL_SPHEROID_A + 2.0*PGL_SPHEROID_B) / 3.0) 10.44 +#define PGL_SCALE (PGL_SPHEROID_A / PGL_DIAMETER) /* semi-major ref. */ 10.45 +#define PGL_FADELIMIT (PGL_DIAMETER * M_PI / 6.0) /* 1/6 circumference */ 10.46 +#define PGL_MAXDIST (PGL_DIAMETER * M_PI / 2.0) /* maximum distance */ 10.47 + 10.48 +/* calculate distance between two points on earth (given in degrees) */ 10.49 +static inline double pgl_distance( 10.50 + double lat1, double lon1, double lat2, double lon2 10.51 +) { 10.52 + float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin; 10.53 + float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t; 10.54 + /* normalize delta longitude (lon2 > 0 && lon1 = 0) */ 10.55 + /* lon1 = 0 (not used anymore) */ 10.56 + lon2 = fabs(lon2-lon1); 10.57 + /* convert to radians (first divide, then multiply) */ 10.58 + lat1 = (lat1 / 180.0) * M_PI; 10.59 + lat2 = (lat2 / 180.0) * M_PI; 10.60 + lon2 = (lon2 / 180.0) * M_PI; 10.61 + /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point 10.62 + operations (lon2 >= lon1 is already ensured in a previous step) */ 10.63 + if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; } 10.64 + /* calculate 3d coordinates on scaled ellipsoid which has an average diameter 10.65 + of 1.0 */ 10.66 + lat1cos = cos(lat1); lat1sin = sin(lat1); 10.67 + lat2cos = cos(lat2); lat2sin = sin(lat2); 10.68 + lon2cos = cos(lon2); lon2sin = sin(lon2); 10.69 + nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin); 10.70 + nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin); 10.71 + x1 = nphi1 * lat1cos; 10.72 + z1 = nphi1 * PGL_SUBEPS2 * lat1sin; 10.73 + x2 = nphi2 * lat2cos * lon2cos; 10.74 + y2 = nphi2 * lat2cos * lon2sin; 10.75 + z2 = nphi2 * PGL_SUBEPS2 * lat2sin; 10.76 + /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */ 10.77 + g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1)); 10.78 + /* convert tunnel distance through scaled ellipsoid to approximated surface 10.79 + distance on original ellipsoid */ 10.80 + if (g > 1.0) g = 1.0; 10.81 + s = PGL_DIAMETER * asin(g); 10.82 + /* return result only if small enough to be precise (less than 1/3 of 10.83 + maximum possible distance) */ 10.84 + if (s <= PGL_FADELIMIT) return s; 10.85 + /* determine antipodal point of second point (i.e. mirror second point) */ 10.86 + lat2 = -lat2; lon2 = lon2 - M_PI; 10.87 + lat2cos = cos(lat2); lat2sin = sin(lat2); 10.88 + lon2cos = cos(lon2); lon2sin = sin(lon2); 10.89 + /* calculate 3d coordinates of antipodal point on scaled ellipsoid */ 10.90 + nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin); 10.91 + x2 = nphi2 * lat2cos * lon2cos; 10.92 + y2 = nphi2 * lat2cos * lon2sin; 10.93 + z2 = nphi2 * PGL_SUBEPS2 * lat2sin; 10.94 + /* calculate tunnel distance to antipodal point through scaled ellipsoid */ 10.95 + g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1)); 10.96 + /* convert tunnel distance to antipodal point through scaled ellipsoid to 10.97 + approximated surface distance to antipodal point on original ellipsoid */ 10.98 + if (g > 1.0) g = 1.0; 10.99 + t = PGL_DIAMETER * asin(g); 10.100 + /* surface distance between original points can now be approximated by 10.101 + substracting antipodal distance from maximum possible distance; 10.102 + return result only if small enough (less than 1/3 of maximum possible 10.103 + distance) */ 10.104 + if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t; 10.105 + /* otherwise crossfade direct and antipodal result to ensure monotonicity */ 10.106 + return ( 10.107 + (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) / 10.108 + (s + t - 2*PGL_FADELIMIT) 10.109 + ); 10.110 +} 10.111 + 10.112 +/* finite distance that can not be reached on earth */ 10.113 +#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST) 10.114 + 10.115 + 10.116 +/*--------------------------------* 10.117 + * simple geographic data types * 10.118 + *--------------------------------*/ 10.119 + 10.120 +/* point on earth given by latitude and longitude in degrees */ 10.121 +/* (type "epoint" in SQL) */ 10.122 +typedef struct { 10.123 + double lat; /* between -90 and 90 (both inclusive) */ 10.124 + double lon; /* between -180 and 180 (both inclusive) */ 10.125 +} pgl_point; 10.126 + 10.127 +/* box delimited by two parallels and two meridians (all in degrees) */ 10.128 +/* (type "ebox" in SQL) */ 10.129 +typedef struct { 10.130 + double lat_min; /* between -90 and 90 (both inclusive) */ 10.131 + double lat_max; /* between -90 and 90 (both inclusive) */ 10.132 + double lon_min; /* between -180 and 180 (both inclusive) */ 10.133 + double lon_max; /* between -180 and 180 (both inclusive) */ 10.134 + /* if lat_min > lat_max, then box is empty */ 10.135 + /* if lon_min > lon_max, then 180th meridian is crossed */ 10.136 +} pgl_box; 10.137 + 10.138 +/* circle on earth surface (for radial searches with fixed radius) */ 10.139 +/* (type "ecircle" in SQL) */ 10.140 +typedef struct { 10.141 + pgl_point center; 10.142 + double radius; /* positive (including +0 but excluding -0), or -INFINITY */ 10.143 + /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point), 10.144 + zero radius (0) denotes a single point, 10.145 + a finite radius (0 < radius < INFINITY) denotes a filled circle, and 10.146 + a radius of INFINITY is valid and means complete coverage of earth. */ 10.147 +} pgl_circle; 10.148 + 10.149 + 10.150 +/*----------------------------------* 10.151 + * geographic "cluster" data type * 10.152 + *----------------------------------*/ 10.153 + 10.154 +/* A cluster is a collection of points, paths, outlines, and polygons. If two 10.155 + polygons in a cluster overlap, the area covered by both polygons does not 10.156 + belong to the cluster. This way, a cluster can be used to describe complex 10.157 + shapes like polygons with holes. Outlines are non-filled polygons. Paths are 10.158 + open by default (i.e. the last point in the list is not connected with the 10.159 + first point in the list). Note that each outline or polygon in a cluster 10.160 + must cover a longitude range of less than 180 degrees to avoid ambiguities. 10.161 + Areas which are larger may be split into multiple polygons. */ 10.162 + 10.163 +/* maximum number of points in a cluster */ 10.164 +/* (limited to avoid integer overflows, e.g. when allocating memory) */ 10.165 +#define PGL_CLUSTER_MAXPOINTS 16777216 10.166 + 10.167 +/* types of cluster entries */ 10.168 +#define PGL_ENTRY_POINT 1 /* a point */ 10.169 +#define PGL_ENTRY_PATH 2 /* a path from first point to last point */ 10.170 +#define PGL_ENTRY_OUTLINE 3 /* a non-filled polygon with given vertices */ 10.171 +#define PGL_ENTRY_POLYGON 4 /* a filled polygon with given vertices */ 10.172 + 10.173 +/* Entries of a cluster are described by two different structs: pgl_newentry 10.174 + and pgl_entry. The first is used only during construction of a cluster, the 10.175 + second is used in all other cases (e.g. when reading clusters from the 10.176 + database, performing operations, etc). */ 10.177 + 10.178 +/* entry for new geographic cluster during construction of that cluster */ 10.179 +typedef struct { 10.180 + int32_t entrytype; 10.181 + int32_t npoints; 10.182 + pgl_point *points; /* pointer to an array of points (pgl_point) */ 10.183 +} pgl_newentry; 10.184 + 10.185 +/* entry of geographic cluster */ 10.186 +typedef struct { 10.187 + int32_t entrytype; /* type of entry: point, path, outline, polygon */ 10.188 + int32_t npoints; /* number of stored points (set to 1 for point entry) */ 10.189 + int32_t offset; /* offset of pgl_point array from cluster base address */ 10.190 + /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */ 10.191 +} pgl_entry; 10.192 + 10.193 +/* geographic cluster which is a collection of points, (open) paths, polygons, 10.194 + and outlines (non-filled polygons) */ 10.195 +typedef struct { 10.196 + char header[VARHDRSZ]; /* PostgreSQL header for variable size data types */ 10.197 + int32_t nentries; /* number of stored points */ 10.198 + pgl_circle bounding; /* bounding circle */ 10.199 + /* Note: bounding circle ensures alignment of pgl_cluster for points */ 10.200 + pgl_entry entries[FLEXIBLE_ARRAY_MEMBER]; /* var-length data */ 10.201 +} pgl_cluster; 10.202 + 10.203 +/* macro to determine memory alignment of points */ 10.204 +/* (needed to store pgl_point array after entries in pgl_cluster) */ 10.205 +typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment; 10.206 +#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned) 10.207 + 10.208 +/* macro to extract a pointer to the array of points of a cluster entry */ 10.209 +#define PGL_ENTRY_POINTS(cluster, idx) \ 10.210 + ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset)) 10.211 + 10.212 +/* convert pgl_newentry array to pgl_cluster */ 10.213 +static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) { 10.214 + int i; /* index of current entry */ 10.215 + int npoints = 0; /* number of points in whole cluster */ 10.216 + int entry_npoints; /* number of points in current entry */ 10.217 + int points_offset = PGL_POINT_ALIGNMENT * ( 10.218 + ( offsetof(pgl_cluster, entries) + 10.219 + nentries * sizeof(pgl_entry) + 10.220 + PGL_POINT_ALIGNMENT - 1 10.221 + ) / PGL_POINT_ALIGNMENT 10.222 + ); /* offset of pgl_point array from base address (considering alignment) */ 10.223 + pgl_cluster *cluster; /* new cluster to be returned */ 10.224 + /* determine total number of points */ 10.225 + for (i=0; i<nentries; i++) npoints += entries[i].npoints; 10.226 + /* allocate memory for cluster (including entries and points) */ 10.227 + cluster = palloc(points_offset + npoints * sizeof(pgl_point)); 10.228 + /* re-count total number of points to determine offset for each entry */ 10.229 + npoints = 0; 10.230 + /* copy entries and points */ 10.231 + for (i=0; i<nentries; i++) { 10.232 + /* determine number of points in entry */ 10.233 + entry_npoints = entries[i].npoints; 10.234 + /* copy entry */ 10.235 + cluster->entries[i].entrytype = entries[i].entrytype; 10.236 + cluster->entries[i].npoints = entry_npoints; 10.237 + /* calculate offset (in bytes) of pgl_point array */ 10.238 + cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point); 10.239 + /* copy points */ 10.240 + memcpy( 10.241 + PGL_ENTRY_POINTS(cluster, i), 10.242 + entries[i].points, 10.243 + entry_npoints * sizeof(pgl_point) 10.244 + ); 10.245 + /* update total number of points processed */ 10.246 + npoints += entry_npoints; 10.247 + } 10.248 + /* set number of entries in cluster */ 10.249 + cluster->nentries = nentries; 10.250 + /* set PostgreSQL header for variable sized data */ 10.251 + SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point)); 10.252 + /* return newly created cluster */ 10.253 + return cluster; 10.254 +} 10.255 + 10.256 + 10.257 +/*----------------------------------------* 10.258 + * C functions on geographic data types * 10.259 + *----------------------------------------*/ 10.260 + 10.261 +/* round latitude or longitude to 12 digits after decimal point */ 10.262 +static inline double pgl_round(double val) { 10.263 + return round(val * 1e12) / 1e12; 10.264 +} 10.265 + 10.266 +/* compare two points */ 10.267 +/* (equality when same point on earth is described, otherwise an arbitrary 10.268 + linear order) */ 10.269 +static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) { 10.270 + double lon1, lon2; /* modified longitudes for special cases */ 10.271 + /* use latitude as first ordering criterion */ 10.272 + if (point1->lat < point2->lat) return -1; 10.273 + if (point1->lat > point2->lat) return 1; 10.274 + /* determine modified longitudes (considering special case of poles and 10.275 + 180th meridian which can be described as W180 or E180) */ 10.276 + if (point1->lat == -90 || point1->lat == 90) lon1 = 0; 10.277 + else if (point1->lon == 180) lon1 = -180; 10.278 + else lon1 = point1->lon; 10.279 + if (point2->lat == -90 || point2->lat == 90) lon2 = 0; 10.280 + else if (point2->lon == 180) lon2 = -180; 10.281 + else lon2 = point2->lon; 10.282 + /* use (modified) longitude as secondary ordering criterion */ 10.283 + if (lon1 < lon2) return -1; 10.284 + if (lon1 > lon2) return 1; 10.285 + /* no difference found, points are equal */ 10.286 + return 0; 10.287 +} 10.288 + 10.289 +/* compare two boxes */ 10.290 +/* (equality when same box on earth is described, otherwise an arbitrary linear 10.291 + order) */ 10.292 +static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) { 10.293 + /* two empty boxes are equal, and an empty box is always considered "less 10.294 + than" a non-empty box */ 10.295 + if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1; 10.296 + if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0; 10.297 + if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1; 10.298 + /* use southern border as first ordering criterion */ 10.299 + if (box1->lat_min < box2->lat_min) return -1; 10.300 + if (box1->lat_min > box2->lat_min) return 1; 10.301 + /* use northern border as second ordering criterion */ 10.302 + if (box1->lat_max < box2->lat_max) return -1; 10.303 + if (box1->lat_max > box2->lat_max) return 1; 10.304 + /* use western border as third ordering criterion */ 10.305 + if (box1->lon_min < box2->lon_min) return -1; 10.306 + if (box1->lon_min > box2->lon_min) return 1; 10.307 + /* use eastern border as fourth ordering criterion */ 10.308 + if (box1->lon_max < box2->lon_max) return -1; 10.309 + if (box1->lon_max > box2->lon_max) return 1; 10.310 + /* no difference found, boxes are equal */ 10.311 + return 0; 10.312 +} 10.313 + 10.314 +/* compare two circles */ 10.315 +/* (equality when same circle on earth is described, otherwise an arbitrary 10.316 + linear order) */ 10.317 +static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) { 10.318 + /* two circles with same infinite radius (positive or negative infinity) are 10.319 + considered equal independently of center point */ 10.320 + if ( 10.321 + !isfinite(circle1->radius) && !isfinite(circle2->radius) && 10.322 + circle1->radius == circle2->radius 10.323 + ) return 0; 10.324 + /* use radius as first ordering criterion */ 10.325 + if (circle1->radius < circle2->radius) return -1; 10.326 + if (circle1->radius > circle2->radius) return 1; 10.327 + /* use center point as secondary ordering criterion */ 10.328 + return pgl_point_cmp(&(circle1->center), &(circle2->center)); 10.329 +} 10.330 + 10.331 +/* set box to empty box*/ 10.332 +static void pgl_box_set_empty(pgl_box *box) { 10.333 + box->lat_min = INFINITY; 10.334 + box->lat_max = -INFINITY; 10.335 + box->lon_min = 0; 10.336 + box->lon_max = 0; 10.337 +} 10.338 + 10.339 +/* check if point is inside a box */ 10.340 +static bool pgl_point_in_box(pgl_point *point, pgl_box *box) { 10.341 + return ( 10.342 + point->lat >= box->lat_min && point->lat <= box->lat_max && ( 10.343 + (box->lon_min > box->lon_max) ? ( 10.344 + /* box crosses 180th meridian */ 10.345 + point->lon >= box->lon_min || point->lon <= box->lon_max 10.346 + ) : ( 10.347 + /* box does not cross the 180th meridian */ 10.348 + point->lon >= box->lon_min && point->lon <= box->lon_max 10.349 + ) 10.350 + ) 10.351 + ); 10.352 +} 10.353 + 10.354 +/* check if two boxes overlap */ 10.355 +static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) { 10.356 + return ( 10.357 + box2->lat_max >= box2->lat_min && /* ensure box2 is not empty */ 10.358 + ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) && 10.359 + ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && ( 10.360 + ( 10.361 + /* check if one and only one box crosses the 180th meridian */ 10.362 + ((box1->lon_min > box1->lon_max) ? 1 : 0) ^ 10.363 + ((box2->lon_min > box2->lon_max) ? 1 : 0) 10.364 + ) ? ( 10.365 + /* exactly one box crosses the 180th meridian */ 10.366 + box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min || 10.367 + box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max 10.368 + ) : ( 10.369 + /* no box or both boxes cross the 180th meridian */ 10.370 + ( 10.371 + (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) && 10.372 + (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max) 10.373 + ) || 10.374 + /* handle W180 == E180 */ 10.375 + ( box1->lon_min == -180 && box2->lon_max == 180 ) || 10.376 + ( box2->lon_min == -180 && box1->lon_max == 180 ) 10.377 + ) 10.378 + ) 10.379 + ); 10.380 +} 10.381 + 10.382 +/* check unambiguousness of east/west orientation of cluster entries and set 10.383 + bounding circle of cluster */ 10.384 +static bool pgl_finalize_cluster(pgl_cluster *cluster) { 10.385 + int i, j; /* i: index of entry, j: index of point in entry */ 10.386 + int npoints; /* number of points in entry */ 10.387 + int total_npoints = 0; /* total number of points in cluster */ 10.388 + pgl_point *points; /* points in entry */ 10.389 + int lon_dir; /* first point of entry west (-1) or east (+1) */ 10.390 + double lon_break = 0; /* antipodal longitude of first point in entry */ 10.391 + double lon_min, lon_max; /* covered longitude range of entry */ 10.392 + double value; /* temporary variable */ 10.393 + /* reset bounding circle center to empty circle at 0/0 coordinates */ 10.394 + cluster->bounding.center.lat = 0; 10.395 + cluster->bounding.center.lon = 0; 10.396 + cluster->bounding.radius = -INFINITY; 10.397 + /* if cluster is not empty */ 10.398 + if (cluster->nentries != 0) { 10.399 + /* iterate over all cluster entries and ensure they each cover a longitude 10.400 + range less than 180 degrees */ 10.401 + for (i=0; i<cluster->nentries; i++) { 10.402 + /* get properties of entry */ 10.403 + npoints = cluster->entries[i].npoints; 10.404 + points = PGL_ENTRY_POINTS(cluster, i); 10.405 + /* get longitude of first point of entry */ 10.406 + value = points[0].lon; 10.407 + /* initialize lon_min and lon_max with longitude of first point */ 10.408 + lon_min = value; 10.409 + lon_max = value; 10.410 + /* determine east/west orientation of first point and calculate antipodal 10.411 + longitude (Note: rounding required here) */ 10.412 + if (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); } 10.413 + else if (value > 0) { lon_dir = 1; lon_break = pgl_round(value - 180); } 10.414 + else lon_dir = 0; 10.415 + /* iterate over all other points in entry */ 10.416 + for (j=1; j<npoints; j++) { 10.417 + /* consider longitude wrap-around */ 10.418 + value = points[j].lon; 10.419 + if (lon_dir<0 && value>lon_break) value = pgl_round(value - 360); 10.420 + else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360); 10.421 + /* update lon_min and lon_max */ 10.422 + if (value < lon_min) lon_min = value; 10.423 + else if (value > lon_max) lon_max = value; 10.424 + /* return false if 180 degrees or more are covered */ 10.425 + if (lon_max - lon_min >= 180) return false; 10.426 + } 10.427 + } 10.428 + /* iterate over all points of all entries and calculate arbitrary center 10.429 + point for bounding circle (best if center point minimizes the radius, 10.430 + but some error is allowed here) */ 10.431 + for (i=0; i<cluster->nentries; i++) { 10.432 + /* get properties of entry */ 10.433 + npoints = cluster->entries[i].npoints; 10.434 + points = PGL_ENTRY_POINTS(cluster, i); 10.435 + /* check if first entry */ 10.436 + if (i==0) { 10.437 + /* get longitude of first point of first entry in whole cluster */ 10.438 + value = points[0].lon; 10.439 + /* initialize lon_min and lon_max with longitude of first point of 10.440 + first entry in whole cluster (used to determine if whole cluster 10.441 + covers a longitude range of 180 degrees or more) */ 10.442 + lon_min = value; 10.443 + lon_max = value; 10.444 + /* determine east/west orientation of first point and calculate 10.445 + antipodal longitude (Note: rounding not necessary here) */ 10.446 + if (value < 0) { lon_dir = -1; lon_break = value + 180; } 10.447 + else if (value > 0) { lon_dir = 1; lon_break = value - 180; } 10.448 + else lon_dir = 0; 10.449 + } 10.450 + /* iterate over all points in entry */ 10.451 + for (j=0; j<npoints; j++) { 10.452 + /* longitude wrap-around (Note: rounding not necessary here) */ 10.453 + value = points[j].lon; 10.454 + if (lon_dir < 0 && value > lon_break) value -= 360; 10.455 + else if (lon_dir > 0 && value < lon_break) value += 360; 10.456 + if (value < lon_min) lon_min = value; 10.457 + else if (value > lon_max) lon_max = value; 10.458 + /* set bounding circle to cover whole earth if more than 180 degrees 10.459 + are covered */ 10.460 + if (lon_max - lon_min >= 180) { 10.461 + cluster->bounding.center.lat = 0; 10.462 + cluster->bounding.center.lon = 0; 10.463 + cluster->bounding.radius = INFINITY; 10.464 + return true; 10.465 + } 10.466 + /* add point to bounding circle center (for average calculation) */ 10.467 + cluster->bounding.center.lat += points[j].lat; 10.468 + cluster->bounding.center.lon += value; 10.469 + } 10.470 + /* count total number of points */ 10.471 + total_npoints += npoints; 10.472 + } 10.473 + /* determine average latitude and longitude of cluster */ 10.474 + cluster->bounding.center.lat /= total_npoints; 10.475 + cluster->bounding.center.lon /= total_npoints; 10.476 + /* normalize longitude of center of cluster bounding circle */ 10.477 + if (cluster->bounding.center.lon < -180) { 10.478 + cluster->bounding.center.lon += 360; 10.479 + } 10.480 + else if (cluster->bounding.center.lon > 180) { 10.481 + cluster->bounding.center.lon -= 360; 10.482 + } 10.483 + /* round bounding circle center (useful if it is used by other functions) */ 10.484 + cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat); 10.485 + cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon); 10.486 + /* calculate radius of bounding circle */ 10.487 + for (i=0; i<cluster->nentries; i++) { 10.488 + npoints = cluster->entries[i].npoints; 10.489 + points = PGL_ENTRY_POINTS(cluster, i); 10.490 + for (j=0; j<npoints; j++) { 10.491 + value = pgl_distance( 10.492 + cluster->bounding.center.lat, cluster->bounding.center.lon, 10.493 + points[j].lat, points[j].lon 10.494 + ); 10.495 + if (value > cluster->bounding.radius) cluster->bounding.radius = value; 10.496 + } 10.497 + } 10.498 + } 10.499 + /* return true (east/west orientation is unambiguous) */ 10.500 + return true; 10.501 +} 10.502 + 10.503 +/* check if point is inside cluster */ 10.504 +static bool pgl_point_in_cluster(pgl_point *point, pgl_cluster *cluster) { 10.505 + int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 10.506 + int entrytype; /* type of entry */ 10.507 + int npoints; /* number of points in entry */ 10.508 + pgl_point *points; /* array of points in entry */ 10.509 + int lon_dir = 0; /* first vertex west (-1) or east (+1) */ 10.510 + double lon_break = 0; /* antipodal longitude of first vertex */ 10.511 + double lat0 = point->lat; /* latitude of point */ 10.512 + double lon0; /* (adjusted) longitude of point */ 10.513 + double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ 10.514 + double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ 10.515 + double lon; /* longitude of intersection */ 10.516 + int counter = 0; /* counter for intersections east of point */ 10.517 + /* points outside bounding circle are always assumed to be non-overlapping */ 10.518 + /* (necessary for consistent table and index scans) */ 10.519 + if ( 10.520 + pgl_distance( 10.521 + point->lat, point->lon, 10.522 + cluster->bounding.center.lat, cluster->bounding.center.lon 10.523 + ) > cluster->bounding.radius 10.524 + ) return false; 10.525 + /* iterate over all entries */ 10.526 + for (i=0; i<cluster->nentries; i++) { 10.527 + /* get properties of entry */ 10.528 + entrytype = cluster->entries[i].entrytype; 10.529 + npoints = cluster->entries[i].npoints; 10.530 + points = PGL_ENTRY_POINTS(cluster, i); 10.531 + /* determine east/west orientation of first point of entry and calculate 10.532 + antipodal longitude */ 10.533 + lon_break = points[0].lon; 10.534 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 10.535 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 10.536 + else lon_dir = 0; 10.537 + /* get longitude of point */ 10.538 + lon0 = point->lon; 10.539 + /* consider longitude wrap-around for point */ 10.540 + if (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360); 10.541 + else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360); 10.542 + /* iterate over all edges and vertices */ 10.543 + for (j=0; j<npoints; j++) { 10.544 + /* return true if point is on vertex of polygon */ 10.545 + if (pgl_point_cmp(point, &(points[j])) == 0) return true; 10.546 + /* calculate index of next vertex */ 10.547 + k = (j+1) % npoints; 10.548 + /* skip last edge unless entry is (closed) outline or polygon */ 10.549 + if ( 10.550 + k == 0 && 10.551 + entrytype != PGL_ENTRY_OUTLINE && 10.552 + entrytype != PGL_ENTRY_POLYGON 10.553 + ) continue; 10.554 + /* get latitude and longitude values of edge */ 10.555 + lat1 = points[j].lat; 10.556 + lat2 = points[k].lat; 10.557 + lon1 = points[j].lon; 10.558 + lon2 = points[k].lon; 10.559 + /* consider longitude wrap-around for edge */ 10.560 + if (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360); 10.561 + else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360); 10.562 + if (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360); 10.563 + else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360); 10.564 + /* return true if point is on horizontal (west to east) edge of polygon */ 10.565 + if ( 10.566 + lat0 == lat1 && lat0 == lat2 && 10.567 + ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) ) 10.568 + ) return true; 10.569 + /* check if edge crosses east/west line of point */ 10.570 + if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) { 10.571 + /* calculate longitude of intersection */ 10.572 + lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1); 10.573 + /* return true if intersection goes (approximately) through point */ 10.574 + if (pgl_round(lon) == lon0) return true; 10.575 + /* count intersection if east of point and entry is polygon*/ 10.576 + if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++; 10.577 + } 10.578 + } 10.579 + } 10.580 + /* return true if number of intersections is odd */ 10.581 + return counter & 1; 10.582 +} 10.583 + 10.584 +/* calculate (approximate) distance between point and cluster */ 10.585 +static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) { 10.586 + int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 10.587 + int entrytype; /* type of entry */ 10.588 + int npoints; /* number of points in entry */ 10.589 + pgl_point *points; /* array of points in entry */ 10.590 + int lon_dir = 0; /* first vertex west (-1) or east (+1) */ 10.591 + double lon_break = 0; /* antipodal longitude of first vertex */ 10.592 + double lon_min = 0; /* minimum (adjusted) longitude of entry vertices */ 10.593 + double lon_max = 0; /* maximum (adjusted) longitude of entry vertices */ 10.594 + double lat0 = point->lat; /* latitude of point */ 10.595 + double lon0; /* (adjusted) longitude of point */ 10.596 + double lat1, lon1; /* latitude and (adjusted) longitude of vertex */ 10.597 + double lat2, lon2; /* latitude and (adjusted) longitude of next vertex */ 10.598 + double s; /* scalar for vector calculations */ 10.599 + double dist; /* distance calculated in one step */ 10.600 + double min_dist = INFINITY; /* minimum distance */ 10.601 + /* distance is zero if point is contained in cluster */ 10.602 + if (pgl_point_in_cluster(point, cluster)) return 0; 10.603 + /* iterate over all entries */ 10.604 + for (i=0; i<cluster->nentries; i++) { 10.605 + /* get properties of entry */ 10.606 + entrytype = cluster->entries[i].entrytype; 10.607 + npoints = cluster->entries[i].npoints; 10.608 + points = PGL_ENTRY_POINTS(cluster, i); 10.609 + /* determine east/west orientation of first point of entry and calculate 10.610 + antipodal longitude */ 10.611 + lon_break = points[0].lon; 10.612 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 10.613 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 10.614 + else lon_dir = 0; 10.615 + /* determine covered longitude range */ 10.616 + for (j=0; j<npoints; j++) { 10.617 + /* get longitude of vertex */ 10.618 + lon1 = points[j].lon; 10.619 + /* adjust longitude to fix potential wrap-around */ 10.620 + if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360; 10.621 + else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; 10.622 + /* update minimum and maximum longitude of polygon */ 10.623 + if (j == 0 || lon1 < lon_min) lon_min = lon1; 10.624 + if (j == 0 || lon1 > lon_max) lon_max = lon1; 10.625 + } 10.626 + /* adjust longitude wrap-around according to full longitude range */ 10.627 + lon_break = (lon_max + lon_min) / 2; 10.628 + if (lon_break < 0) { lon_dir = -1; lon_break += 180; } 10.629 + else if (lon_break > 0) { lon_dir = 1; lon_break -= 180; } 10.630 + /* get longitude of point */ 10.631 + lon0 = point->lon; 10.632 + /* consider longitude wrap-around for point */ 10.633 + if (lon_dir < 0 && lon0 > lon_break) lon0 -= 360; 10.634 + else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360; 10.635 + /* iterate over all edges and vertices */ 10.636 + for (j=0; j<npoints; j++) { 10.637 + /* get latitude and longitude values of current point */ 10.638 + lat1 = points[j].lat; 10.639 + lon1 = points[j].lon; 10.640 + /* consider longitude wrap-around for current point */ 10.641 + if (lon_dir < 0 && lon1 > lon_break) lon1 -= 360; 10.642 + else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360; 10.643 + /* calculate distance to vertex */ 10.644 + dist = pgl_distance(lat0, lon0, lat1, lon1); 10.645 + /* store calculated distance if smallest */ 10.646 + if (dist < min_dist) min_dist = dist; 10.647 + /* calculate index of next vertex */ 10.648 + k = (j+1) % npoints; 10.649 + /* skip last edge unless entry is (closed) outline or polygon */ 10.650 + if ( 10.651 + k == 0 && 10.652 + entrytype != PGL_ENTRY_OUTLINE && 10.653 + entrytype != PGL_ENTRY_POLYGON 10.654 + ) continue; 10.655 + /* get latitude and longitude values of next point */ 10.656 + lat2 = points[k].lat; 10.657 + lon2 = points[k].lon; 10.658 + /* consider longitude wrap-around for next point */ 10.659 + if (lon_dir < 0 && lon2 > lon_break) lon2 -= 360; 10.660 + else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360; 10.661 + /* go to next vertex and edge if edge is degenerated */ 10.662 + if (lat1 == lat2 && lon1 == lon2) continue; 10.663 + /* otherwise test if point can be projected onto edge of polygon */ 10.664 + s = ( 10.665 + ((lat0-lat1) * (lat2-lat1) + (lon0-lon1) * (lon2-lon1)) / 10.666 + ((lat2-lat1) * (lat2-lat1) + (lon2-lon1) * (lon2-lon1)) 10.667 + ); 10.668 + /* go to next vertex and edge if point cannot be projected */ 10.669 + if (!(s > 0 && s < 1)) continue; 10.670 + /* calculate distance from original point to projected point */ 10.671 + dist = pgl_distance( 10.672 + lat0, lon0, 10.673 + lat1 + s * (lat2-lat1), 10.674 + lon1 + s * (lon2-lon1) 10.675 + ); 10.676 + /* store calculated distance if smallest */ 10.677 + if (dist < min_dist) min_dist = dist; 10.678 + } 10.679 + } 10.680 + /* return minimum distance */ 10.681 + return min_dist; 10.682 +} 10.683 + 10.684 +/* estimator function for distance between box and point */ 10.685 +/* allowed to return smaller values than actually correct */ 10.686 +static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) { 10.687 + double dlon; /* longitude range of box (delta longitude) */ 10.688 + double h; /* half of distance along meridian */ 10.689 + double d; /* distance between both southern or both northern points */ 10.690 + double cur_dist; /* calculated distance */ 10.691 + double min_dist; /* minimum distance calculated */ 10.692 + /* return infinity if bounding box is empty */ 10.693 + if (box->lat_min > box->lat_max) return INFINITY; 10.694 + /* return zero if point is inside bounding box */ 10.695 + if (pgl_point_in_box(point, box)) return 0; 10.696 + /* calculate delta longitude */ 10.697 + dlon = box->lon_max - box->lon_min; 10.698 + if (dlon < 0) dlon += 360; /* 180th meridian crossed */ 10.699 + /* if delta longitude is greater than 180 degrees, perform safe fall-back */ 10.700 + if (dlon > 180) return 0; 10.701 + /* calculate half of distance along meridian */ 10.702 + h = pgl_distance(box->lat_min, 0, box->lat_max, 0) / 2; 10.703 + /* calculate full distance between southern points */ 10.704 + d = pgl_distance(box->lat_min, 0, box->lat_min, dlon); 10.705 + /* calculate maximum of full distance and half distance */ 10.706 + if (h > d) d = h; 10.707 + /* calculate distance from point to first southern vertex and substract 10.708 + maximum error */ 10.709 + min_dist = pgl_distance( 10.710 + point->lat, point->lon, box->lat_min, box->lon_min 10.711 + ) - d; 10.712 + /* return zero if estimated distance is smaller than zero */ 10.713 + if (min_dist <= 0) return 0; 10.714 + /* repeat procedure with second southern vertex */ 10.715 + cur_dist = pgl_distance( 10.716 + point->lat, point->lon, box->lat_min, box->lon_max 10.717 + ) - d; 10.718 + if (cur_dist <= 0) return 0; 10.719 + if (cur_dist < min_dist) min_dist = cur_dist; 10.720 + /* calculate full distance between northern points */ 10.721 + d = pgl_distance(box->lat_max, 0, box->lat_max, dlon); 10.722 + /* calculate maximum of full distance and half distance */ 10.723 + if (h > d) d = h; 10.724 + /* repeat procedure with northern vertices */ 10.725 + cur_dist = pgl_distance( 10.726 + point->lat, point->lon, box->lat_max, box->lon_max 10.727 + ) - d; 10.728 + if (cur_dist <= 0) return 0; 10.729 + if (cur_dist < min_dist) min_dist = cur_dist; 10.730 + cur_dist = pgl_distance( 10.731 + point->lat, point->lon, box->lat_max, box->lon_min 10.732 + ) - d; 10.733 + if (cur_dist <= 0) return 0; 10.734 + if (cur_dist < min_dist) min_dist = cur_dist; 10.735 + /* return smallest value (unless already returned zero) */ 10.736 + return min_dist; 10.737 +} 10.738 + 10.739 + 10.740 +/*----------------------------* 10.741 + * fractal geographic index * 10.742 + *----------------------------*/ 10.743 + 10.744 +/* number of bytes used for geographic (center) position in keys */ 10.745 +#define PGL_KEY_LATLON_BYTELEN 7 10.746 + 10.747 +/* maximum reference value for logarithmic size of geographic objects */ 10.748 +#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0) /* can be tweaked */ 10.749 + 10.750 +/* safety margin to avoid floating point errors in distance estimation */ 10.751 +#define PGL_FPE_SAFETY (1.0+1e-14) /* slightly greater than 1.0 */ 10.752 + 10.753 +/* pointer to index key (either pgl_pointkey or pgl_areakey) */ 10.754 +typedef unsigned char *pgl_keyptr; 10.755 + 10.756 +/* index key for points (objects with zero area) on the spheroid */ 10.757 +/* bit 0..55: interspersed bits of latitude and longitude, 10.758 + bit 56..57: always zero, 10.759 + bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */ 10.760 +typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1]; 10.761 + 10.762 +/* index key for geographic objects on spheroid with area greater than zero */ 10.763 +/* bit 0..55: interspersed bits of latitude and longitude of center point, 10.764 + bit 56: always set to 1, 10.765 + bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.), 10.766 + bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to 10.767 + PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth 10.768 + = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL 10.769 + (with interspersed bits = 0 and node depth = 0) for keys which 10.770 + cover both empty and non-empty objects */ 10.771 + 10.772 +typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2]; 10.773 + 10.774 +/* helper macros for reading/writing index keys */ 10.775 +#define PGL_KEY_NODEDEPTH_OFFSET PGL_KEY_LATLON_BYTELEN 10.776 +#define PGL_KEY_OBJSIZE_OFFSET (PGL_KEY_NODEDEPTH_OFFSET+1) 10.777 +#define PGL_POINTKEY_MAXDEPTH (PGL_KEY_LATLON_BYTELEN*8) 10.778 +#define PGL_AREAKEY_MAXDEPTH (2*PGL_POINTKEY_MAXDEPTH+1) 10.779 +#define PGL_AREAKEY_MAXOBJSIZE (PGL_POINTKEY_MAXDEPTH+1) 10.780 +#define PGL_AREAKEY_TYPEMASK 0x80 10.781 +#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8))) 10.782 +#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \ 10.783 + ( PGL_KEY_LATLONBIT(key1, n) ^ \ 10.784 + PGL_KEY_LATLONBIT(key2, n) ) 10.785 +#define PGL_KEY_IS_AREAKEY(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ 10.786 + PGL_AREAKEY_TYPEMASK) 10.787 +#define PGL_KEY_NODEDEPTH(key) ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \ 10.788 + (PGL_AREAKEY_TYPEMASK-1)) 10.789 +#define PGL_KEY_OBJSIZE(key) ((key)[PGL_KEY_OBJSIZE_OFFSET]) 10.790 +#define PGL_KEY_OBJSIZE_EMPTY 126 10.791 +#define PGL_KEY_OBJSIZE_UNIVERSAL 127 10.792 +#define PGL_KEY_IS_EMPTY(key) ( PGL_KEY_IS_AREAKEY(key) && \ 10.793 + (key)[PGL_KEY_OBJSIZE_OFFSET] == \ 10.794 + PGL_KEY_OBJSIZE_EMPTY ) 10.795 +#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \ 10.796 + (key)[PGL_KEY_OBJSIZE_OFFSET] == \ 10.797 + PGL_KEY_OBJSIZE_UNIVERSAL ) 10.798 + 10.799 +/* set area key to match empty objects only */ 10.800 +static void pgl_key_set_empty(pgl_keyptr key) { 10.801 + memset(key, 0, sizeof(pgl_areakey)); 10.802 + /* Note: setting node depth to maximum is required for picksplit function */ 10.803 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; 10.804 + key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY; 10.805 +} 10.806 + 10.807 +/* set area key to match any object (including empty objects) */ 10.808 +static void pgl_key_set_universal(pgl_keyptr key) { 10.809 + memset(key, 0, sizeof(pgl_areakey)); 10.810 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK; 10.811 + key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL; 10.812 +} 10.813 + 10.814 +/* convert a point on earth into a max-depth key to be used in index */ 10.815 +static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) { 10.816 + double lat = point->lat; 10.817 + double lon = point->lon; 10.818 + int i; 10.819 + /* clear latitude and longitude bits */ 10.820 + memset(key, 0, PGL_KEY_LATLON_BYTELEN); 10.821 + /* set node depth to maximum and type bit to zero */ 10.822 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH; 10.823 + /* iterate over all latitude/longitude bit pairs */ 10.824 + for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) { 10.825 + /* determine latitude bit */ 10.826 + if (lat >= 0) { 10.827 + key[i/4] |= 0x80 >> (2*(i%4)); 10.828 + lat *= 2; lat -= 90; 10.829 + } else { 10.830 + lat *= 2; lat += 90; 10.831 + } 10.832 + /* determine longitude bit */ 10.833 + if (lon >= 0) { 10.834 + key[i/4] |= 0x80 >> (2*(i%4)+1); 10.835 + lon *= 2; lon -= 180; 10.836 + } else { 10.837 + lon *= 2; lon += 180; 10.838 + } 10.839 + } 10.840 +} 10.841 + 10.842 +/* convert a circle on earth into a max-depth key to be used in an index */ 10.843 +static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) { 10.844 + /* handle special case of empty circle */ 10.845 + if (circle->radius < 0) { 10.846 + pgl_key_set_empty(key); 10.847 + return; 10.848 + } 10.849 + /* perform same action as for point keys */ 10.850 + pgl_point_to_key(&(circle->center), key); 10.851 + /* but overwrite type and node depth to fit area index key */ 10.852 + key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH; 10.853 + /* check if radius is greater than (or equal to) reference size */ 10.854 + /* (treat equal values as greater values for numerical safety) */ 10.855 + if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) { 10.856 + /* if yes, set logarithmic size to zero */ 10.857 + key[PGL_KEY_OBJSIZE_OFFSET] = 0; 10.858 + } else { 10.859 + /* otherwise, determine logarithmic size iteratively */ 10.860 + /* (one step is equivalent to a factor of sqrt(2)) */ 10.861 + double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2; 10.862 + int objsize = 1; 10.863 + while (objsize < PGL_AREAKEY_MAXOBJSIZE) { 10.864 + /* stop when radius is greater than (or equal to) adjusted reference */ 10.865 + /* (treat equal values as greater values for numerical safety) */ 10.866 + if (circle->radius >= reference) break; 10.867 + reference /= M_SQRT2; 10.868 + objsize++; 10.869 + } 10.870 + /* set logarithmic size to determined value */ 10.871 + key[PGL_KEY_OBJSIZE_OFFSET] = objsize; 10.872 + } 10.873 +} 10.874 + 10.875 +/* check if one key is subkey of another key or vice versa */ 10.876 +static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) { 10.877 + int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ 10.878 + /* determine smallest depth */ 10.879 + int depth1 = PGL_KEY_NODEDEPTH(key1); 10.880 + int depth2 = PGL_KEY_NODEDEPTH(key2); 10.881 + int depth = (depth1 < depth2) ? depth1 : depth2; 10.882 + /* check if keys are area keys (assuming that both keys have same type) */ 10.883 + if (PGL_KEY_IS_AREAKEY(key1)) { 10.884 + int j = 0; /* bit offset for logarithmic object size bits */ 10.885 + int k = 0; /* bit offset for latitude and longitude */ 10.886 + /* fetch logarithmic object size information */ 10.887 + int objsize1 = PGL_KEY_OBJSIZE(key1); 10.888 + int objsize2 = PGL_KEY_OBJSIZE(key2); 10.889 + /* handle special cases for empty objects (universal and empty keys) */ 10.890 + if ( 10.891 + objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL || 10.892 + objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL 10.893 + ) return true; 10.894 + if ( 10.895 + objsize1 == PGL_KEY_OBJSIZE_EMPTY || 10.896 + objsize2 == PGL_KEY_OBJSIZE_EMPTY 10.897 + ) return objsize1 == objsize2; 10.898 + /* iterate through key bits */ 10.899 + for (i=0; i<depth; i++) { 10.900 + /* every second bit is a bit describing the object size */ 10.901 + if (i%2 == 0) { 10.902 + /* check if object size bit is different in both keys (objsize1 and 10.903 + objsize2 describe the minimum index when object size bit is set) */ 10.904 + if ( 10.905 + (objsize1 <= j && objsize2 > j) || 10.906 + (objsize2 <= j && objsize1 > j) 10.907 + ) { 10.908 + /* bit differs, therefore keys are in separate branches */ 10.909 + return false; 10.910 + } 10.911 + /* increase bit counter for object size bits */ 10.912 + j++; 10.913 + } 10.914 + /* all other bits describe latitude and longitude */ 10.915 + else { 10.916 + /* check if bit differs in both keys */ 10.917 + if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) { 10.918 + /* bit differs, therefore keys are in separate branches */ 10.919 + return false; 10.920 + } 10.921 + /* increase bit counter for latitude/longitude bits */ 10.922 + k++; 10.923 + } 10.924 + } 10.925 + } 10.926 + /* if not, keys are point keys */ 10.927 + else { 10.928 + /* iterate through key bits */ 10.929 + for (i=0; i<depth; i++) { 10.930 + /* check if bit differs in both keys */ 10.931 + if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) { 10.932 + /* bit differs, therefore keys are in separate branches */ 10.933 + return false; 10.934 + } 10.935 + } 10.936 + } 10.937 + /* return true because keys are in the same branch */ 10.938 + return true; 10.939 +} 10.940 + 10.941 +/* combine two keys into new key which covers both original keys */ 10.942 +/* (result stored in first argument) */ 10.943 +static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) { 10.944 + int i; /* key bit offset (includes both lat/lon and log. obj. size bits) */ 10.945 + /* determine smallest depth */ 10.946 + int depth1 = PGL_KEY_NODEDEPTH(dst); 10.947 + int depth2 = PGL_KEY_NODEDEPTH(src); 10.948 + int depth = (depth1 < depth2) ? depth1 : depth2; 10.949 + /* check if keys are area keys (assuming that both keys have same type) */ 10.950 + if (PGL_KEY_IS_AREAKEY(dst)) { 10.951 + pgl_areakey dstbuf = { 0, }; /* destination buffer (cleared) */ 10.952 + int j = 0; /* bit offset for logarithmic object size bits */ 10.953 + int k = 0; /* bit offset for latitude and longitude */ 10.954 + /* fetch logarithmic object size information */ 10.955 + int objsize1 = PGL_KEY_OBJSIZE(dst); 10.956 + int objsize2 = PGL_KEY_OBJSIZE(src); 10.957 + /* handle special cases for empty objects (universal and empty keys) */ 10.958 + if ( 10.959 + objsize1 > PGL_AREAKEY_MAXOBJSIZE || 10.960 + objsize2 > PGL_AREAKEY_MAXOBJSIZE 10.961 + ) { 10.962 + if ( 10.963 + objsize1 == PGL_KEY_OBJSIZE_EMPTY && 10.964 + objsize2 == PGL_KEY_OBJSIZE_EMPTY 10.965 + ) pgl_key_set_empty(dst); 10.966 + else pgl_key_set_universal(dst); 10.967 + return; 10.968 + } 10.969 + /* iterate through key bits */ 10.970 + for (i=0; i<depth; i++) { 10.971 + /* every second bit is a bit describing the object size */ 10.972 + if (i%2 == 0) { 10.973 + /* increase bit counter for object size bits first */ 10.974 + /* (handy when setting objsize variable) */ 10.975 + j++; 10.976 + /* check if object size bit is set in neither key */ 10.977 + if (objsize1 >= j && objsize2 >= j) { 10.978 + /* set objsize in destination buffer to indicate that size bit is 10.979 + unset in destination buffer at the current bit position */ 10.980 + dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j; 10.981 + } 10.982 + /* break if object size bit is set in one key only */ 10.983 + else if (objsize1 >= j || objsize2 >= j) break; 10.984 + } 10.985 + /* all other bits describe latitude and longitude */ 10.986 + else { 10.987 + /* break if bit differs in both keys */ 10.988 + if (PGL_KEY_LATLONBIT(dst, k)) { 10.989 + if (!PGL_KEY_LATLONBIT(src, k)) break; 10.990 + /* but set bit in destination buffer if bit is set in both keys */ 10.991 + dstbuf[k/8] |= 0x80 >> (k%8); 10.992 + } else if (PGL_KEY_LATLONBIT(src, k)) break; 10.993 + /* increase bit counter for latitude/longitude bits */ 10.994 + k++; 10.995 + } 10.996 + } 10.997 + /* set common node depth and type bit (type bit = 1) */ 10.998 + dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i; 10.999 + /* copy contents of destination buffer to first key */ 10.1000 + memcpy(dst, dstbuf, sizeof(pgl_areakey)); 10.1001 + } 10.1002 + /* if not, keys are point keys */ 10.1003 + else { 10.1004 + pgl_pointkey dstbuf = { 0, }; /* destination buffer (cleared) */ 10.1005 + /* iterate through key bits */ 10.1006 + for (i=0; i<depth; i++) { 10.1007 + /* break if bit differs in both keys */ 10.1008 + if (PGL_KEY_LATLONBIT(dst, i)) { 10.1009 + if (!PGL_KEY_LATLONBIT(src, i)) break; 10.1010 + /* but set bit in destination buffer if bit is set in both keys */ 10.1011 + dstbuf[i/8] |= 0x80 >> (i%8); 10.1012 + } else if (PGL_KEY_LATLONBIT(src, i)) break; 10.1013 + } 10.1014 + /* set common node depth (type bit = 0) */ 10.1015 + dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i; 10.1016 + /* copy contents of destination buffer to first key */ 10.1017 + memcpy(dst, dstbuf, sizeof(pgl_pointkey)); 10.1018 + } 10.1019 +} 10.1020 + 10.1021 +/* determine center(!) boundaries and radius estimation of index key */ 10.1022 +static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) { 10.1023 + int i; 10.1024 + /* determine node depth */ 10.1025 + int depth = PGL_KEY_NODEDEPTH(key); 10.1026 + /* center point of possible result */ 10.1027 + double lat = 0; 10.1028 + double lon = 0; 10.1029 + /* maximum distance of real center point from key center */ 10.1030 + double dlat = 90; 10.1031 + double dlon = 180; 10.1032 + /* maximum radius of contained objects */ 10.1033 + double radius = 0; /* always return zero for point index keys */ 10.1034 + /* check if key is area key */ 10.1035 + if (PGL_KEY_IS_AREAKEY(key)) { 10.1036 + /* get logarithmic object size */ 10.1037 + int objsize = PGL_KEY_OBJSIZE(key); 10.1038 + /* handle special cases for empty objects (universal and empty keys) */ 10.1039 + if (objsize == PGL_KEY_OBJSIZE_EMPTY) { 10.1040 + pgl_box_set_empty(box); 10.1041 + return 0; 10.1042 + } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) { 10.1043 + box->lat_min = -90; 10.1044 + box->lat_max = 90; 10.1045 + box->lon_min = -180; 10.1046 + box->lon_max = 180; 10.1047 + return 0; /* any value >= 0 would do */ 10.1048 + } 10.1049 + /* calculate maximum possible radius of objects covered by the given key */ 10.1050 + if (objsize == 0) radius = INFINITY; 10.1051 + else { 10.1052 + radius = PGL_AREAKEY_REFOBJSIZE; 10.1053 + while (--objsize) radius /= M_SQRT2; 10.1054 + } 10.1055 + /* iterate over latitude and longitude bits in key */ 10.1056 + /* (every second bit is a latitude or longitude bit) */ 10.1057 + for (i=0; i<depth/2; i++) { 10.1058 + /* check if latitude bit */ 10.1059 + if (i%2 == 0) { 10.1060 + /* cut latitude dimension in half */ 10.1061 + dlat /= 2; 10.1062 + /* increase center latitude if bit is 1, otherwise decrease */ 10.1063 + if (PGL_KEY_LATLONBIT(key, i)) lat += dlat; 10.1064 + else lat -= dlat; 10.1065 + } 10.1066 + /* otherwise longitude bit */ 10.1067 + else { 10.1068 + /* cut longitude dimension in half */ 10.1069 + dlon /= 2; 10.1070 + /* increase center longitude if bit is 1, otherwise decrease */ 10.1071 + if (PGL_KEY_LATLONBIT(key, i)) lon += dlon; 10.1072 + else lon -= dlon; 10.1073 + } 10.1074 + } 10.1075 + } 10.1076 + /* if not, keys are point keys */ 10.1077 + else { 10.1078 + /* iterate over all bits in key */ 10.1079 + for (i=0; i<depth; i++) { 10.1080 + /* check if latitude bit */ 10.1081 + if (i%2 == 0) { 10.1082 + /* cut latitude dimension in half */ 10.1083 + dlat /= 2; 10.1084 + /* increase center latitude if bit is 1, otherwise decrease */ 10.1085 + if (PGL_KEY_LATLONBIT(key, i)) lat += dlat; 10.1086 + else lat -= dlat; 10.1087 + } 10.1088 + /* otherwise longitude bit */ 10.1089 + else { 10.1090 + /* cut longitude dimension in half */ 10.1091 + dlon /= 2; 10.1092 + /* increase center longitude if bit is 1, otherwise decrease */ 10.1093 + if (PGL_KEY_LATLONBIT(key, i)) lon += dlon; 10.1094 + else lon -= dlon; 10.1095 + } 10.1096 + } 10.1097 + } 10.1098 + /* calculate boundaries from center point and remaining dlat and dlon */ 10.1099 + /* (return values through pointer to box) */ 10.1100 + box->lat_min = lat - dlat; 10.1101 + box->lat_max = lat + dlat; 10.1102 + box->lon_min = lon - dlon; 10.1103 + box->lon_max = lon + dlon; 10.1104 + /* return radius (as a function return value) */ 10.1105 + return radius; 10.1106 +} 10.1107 + 10.1108 +/* estimator function for distance between point and index key */ 10.1109 +/* allowed to return smaller values than actually correct */ 10.1110 +static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) { 10.1111 + pgl_box box; /* center(!) bounding box of area index key */ 10.1112 + /* calculate center(!) bounding box and maximum radius of objects covered 10.1113 + by area index key (radius is zero for point index keys) */ 10.1114 + double distance = pgl_key_to_box(key, &box); 10.1115 + /* calculate estimated distance between bounding box of center point of 10.1116 + indexed object and point passed as second argument, then substract maximum 10.1117 + radius of objects covered by index key */ 10.1118 + /* (use PGL_FPE_SAFETY factor to cope with minor floating point errors) */ 10.1119 + distance = ( 10.1120 + pgl_estimate_point_box_distance(point, &box) / PGL_FPE_SAFETY - 10.1121 + distance * PGL_FPE_SAFETY 10.1122 + ); 10.1123 + /* truncate negative results to zero */ 10.1124 + if (distance <= 0) distance = 0; 10.1125 + /* return result */ 10.1126 + return distance; 10.1127 +} 10.1128 + 10.1129 + 10.1130 +/*---------------------------------* 10.1131 + * helper functions for text I/O * 10.1132 + *---------------------------------*/ 10.1133 + 10.1134 +#define PGL_NUMBUFLEN 64 /* buffer size for number to string conversion */ 10.1135 + 10.1136 +/* convert floating point number to string (round-trip safe) */ 10.1137 +static void pgl_print_float(char *buf, double flt) { 10.1138 + /* check if number is integral */ 10.1139 + if (trunc(flt) == flt) { 10.1140 + /* for integral floats use maximum precision */ 10.1141 + snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); 10.1142 + } else { 10.1143 + /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */ 10.1144 + snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt); 10.1145 + if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt); 10.1146 + if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt); 10.1147 + } 10.1148 +} 10.1149 + 10.1150 +/* convert latitude floating point number (in degrees) to string */ 10.1151 +static void pgl_print_lat(char *buf, double lat) { 10.1152 + if (signbit(lat)) { 10.1153 + /* treat negative latitudes (including -0) as south */ 10.1154 + snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat); 10.1155 + } else { 10.1156 + /* treat positive latitudes (including +0) as north */ 10.1157 + snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat); 10.1158 + } 10.1159 +} 10.1160 + 10.1161 +/* convert longitude floating point number (in degrees) to string */ 10.1162 +static void pgl_print_lon(char *buf, double lon) { 10.1163 + if (signbit(lon)) { 10.1164 + /* treat negative longitudes (including -0) as west */ 10.1165 + snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon); 10.1166 + } else { 10.1167 + /* treat positive longitudes (including +0) as east */ 10.1168 + snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon); 10.1169 + } 10.1170 +} 10.1171 + 10.1172 +/* bit masks used as return value of pgl_scan() function */ 10.1173 +#define PGL_SCAN_NONE 0 /* no value has been parsed */ 10.1174 +#define PGL_SCAN_LAT (1<<0) /* latitude has been parsed */ 10.1175 +#define PGL_SCAN_LON (1<<1) /* longitude has been parsed */ 10.1176 +#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON) /* bitwise OR of both */ 10.1177 + 10.1178 +/* parse a coordinate (can be latitude or longitude) */ 10.1179 +static int pgl_scan(char **str, double *lat, double *lon) { 10.1180 + double val; 10.1181 + int len; 10.1182 + if ( 10.1183 + sscanf(*str, " N %lf %n", &val, &len) || 10.1184 + sscanf(*str, " n %lf %n", &val, &len) 10.1185 + ) { 10.1186 + *str += len; *lat = val; return PGL_SCAN_LAT; 10.1187 + } 10.1188 + if ( 10.1189 + sscanf(*str, " S %lf %n", &val, &len) || 10.1190 + sscanf(*str, " s %lf %n", &val, &len) 10.1191 + ) { 10.1192 + *str += len; *lat = -val; return PGL_SCAN_LAT; 10.1193 + } 10.1194 + if ( 10.1195 + sscanf(*str, " E %lf %n", &val, &len) || 10.1196 + sscanf(*str, " e %lf %n", &val, &len) 10.1197 + ) { 10.1198 + *str += len; *lon = val; return PGL_SCAN_LON; 10.1199 + } 10.1200 + if ( 10.1201 + sscanf(*str, " W %lf %n", &val, &len) || 10.1202 + sscanf(*str, " w %lf %n", &val, &len) 10.1203 + ) { 10.1204 + *str += len; *lon = -val; return PGL_SCAN_LON; 10.1205 + } 10.1206 + return PGL_SCAN_NONE; 10.1207 +} 10.1208 + 10.1209 + 10.1210 +/*-----------------* 10.1211 + * SQL functions * 10.1212 + *-----------------*/ 10.1213 + 10.1214 +/* Note: These function names use "epoint", "ebox", etc. notation here instead 10.1215 + of "point", "box", etc. in order to distinguish them from any previously 10.1216 + defined functions. */ 10.1217 + 10.1218 +/* function needed for dummy types and/or not implemented features */ 10.1219 +PG_FUNCTION_INFO_V1(pgl_notimpl); 10.1220 +Datum pgl_notimpl(PG_FUNCTION_ARGS) { 10.1221 + ereport(ERROR, (errmsg("not implemented by pgLatLon"))); 10.1222 +} 10.1223 + 10.1224 +/* set point to latitude and longitude (including checks) */ 10.1225 +static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) { 10.1226 + /* reject infinite or NaN values */ 10.1227 + if (!isfinite(lat) || !isfinite(lon)) { 10.1228 + ereport(ERROR, ( 10.1229 + errcode(ERRCODE_DATA_EXCEPTION), 10.1230 + errmsg("epoint requires finite coordinates") 10.1231 + )); 10.1232 + } 10.1233 + /* check latitude bounds */ 10.1234 + if (lat < -90) { 10.1235 + ereport(WARNING, (errmsg("latitude exceeds south pole"))); 10.1236 + lat = -90; 10.1237 + } else if (lat > 90) { 10.1238 + ereport(WARNING, (errmsg("latitude exceeds north pole"))); 10.1239 + lat = 90; 10.1240 + } 10.1241 + /* check longitude bounds */ 10.1242 + if (lon < -180) { 10.1243 + ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized"))); 10.1244 + lon += 360 - trunc(lon / 360) * 360; 10.1245 + } else if (lon > 180) { 10.1246 + ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized"))); 10.1247 + lon -= 360 + trunc(lon / 360) * 360; 10.1248 + } 10.1249 + /* store rounded latitude/longitude values for round-trip safety */ 10.1250 + point->lat = pgl_round(lat); 10.1251 + point->lon = pgl_round(lon); 10.1252 +} 10.1253 + 10.1254 +/* create point ("epoint" in SQL) from latitude and longitude */ 10.1255 +PG_FUNCTION_INFO_V1(pgl_create_epoint); 10.1256 +Datum pgl_create_epoint(PG_FUNCTION_ARGS) { 10.1257 + pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); 10.1258 + pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1)); 10.1259 + PG_RETURN_POINTER(point); 10.1260 +} 10.1261 + 10.1262 +/* parse point ("epoint" in SQL) */ 10.1263 +/* format: '[NS]<float> [EW]<float>' */ 10.1264 +PG_FUNCTION_INFO_V1(pgl_epoint_in); 10.1265 +Datum pgl_epoint_in(PG_FUNCTION_ARGS) { 10.1266 + char *str = PG_GETARG_CSTRING(0); /* input string */ 10.1267 + char *strptr = str; /* current position within string */ 10.1268 + int done = 0; /* bit mask storing if latitude or longitude was read */ 10.1269 + double lat, lon; /* parsed values as double precision floats */ 10.1270 + pgl_point *point; /* return value (to be palloc'ed) */ 10.1271 + /* parse two floats (each latitude or longitude) separated by white-space */ 10.1272 + done |= pgl_scan(&strptr, &lat, &lon); 10.1273 + if (strptr != str && isspace(strptr[-1])) { 10.1274 + done |= pgl_scan(&strptr, &lat, &lon); 10.1275 + } 10.1276 + /* require end of string, and latitude and longitude parsed successfully */ 10.1277 + if (strptr[0] || done != PGL_SCAN_LATLON) { 10.1278 + ereport(ERROR, ( 10.1279 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 10.1280 + errmsg("invalid input syntax for type epoint: \"%s\"", str) 10.1281 + )); 10.1282 + } 10.1283 + /* allocate memory for result */ 10.1284 + point = (pgl_point *)palloc(sizeof(pgl_point)); 10.1285 + /* set latitude and longitude (and perform checks) */ 10.1286 + pgl_epoint_set_latlon(point, lat, lon); 10.1287 + /* return result */ 10.1288 + PG_RETURN_POINTER(point); 10.1289 +} 10.1290 + 10.1291 +/* create box ("ebox" in SQL) that is empty */ 10.1292 +PG_FUNCTION_INFO_V1(pgl_create_empty_ebox); 10.1293 +Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) { 10.1294 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 10.1295 + pgl_box_set_empty(box); 10.1296 + PG_RETURN_POINTER(box); 10.1297 +} 10.1298 + 10.1299 +/* set box to given boundaries (including checks) */ 10.1300 +static void pgl_ebox_set_boundaries( 10.1301 + pgl_box *box, 10.1302 + double lat_min, double lat_max, double lon_min, double lon_max 10.1303 +) { 10.1304 + /* if minimum latitude is greater than maximum latitude, return empty box */ 10.1305 + if (lat_min > lat_max) { 10.1306 + pgl_box_set_empty(box); 10.1307 + return; 10.1308 + } 10.1309 + /* otherwise reject infinite or NaN values */ 10.1310 + if ( 10.1311 + !isfinite(lat_min) || !isfinite(lat_max) || 10.1312 + !isfinite(lon_min) || !isfinite(lon_max) 10.1313 + ) { 10.1314 + ereport(ERROR, ( 10.1315 + errcode(ERRCODE_DATA_EXCEPTION), 10.1316 + errmsg("ebox requires finite coordinates") 10.1317 + )); 10.1318 + } 10.1319 + /* check latitude bounds */ 10.1320 + if (lat_max < -90) { 10.1321 + ereport(WARNING, (errmsg("northern latitude exceeds south pole"))); 10.1322 + lat_max = -90; 10.1323 + } else if (lat_max > 90) { 10.1324 + ereport(WARNING, (errmsg("northern latitude exceeds north pole"))); 10.1325 + lat_max = 90; 10.1326 + } 10.1327 + if (lat_min < -90) { 10.1328 + ereport(WARNING, (errmsg("southern latitude exceeds south pole"))); 10.1329 + lat_min = -90; 10.1330 + } else if (lat_min > 90) { 10.1331 + ereport(WARNING, (errmsg("southern latitude exceeds north pole"))); 10.1332 + lat_min = 90; 10.1333 + } 10.1334 + /* check if all longitudes are included */ 10.1335 + if (lon_max - lon_min >= 360) { 10.1336 + if (lon_max - lon_min > 360) ereport(WARNING, ( 10.1337 + errmsg("longitude coverage greater than 360 degrees") 10.1338 + )); 10.1339 + lon_min = -180; 10.1340 + lon_max = 180; 10.1341 + } else { 10.1342 + /* normalize longitude bounds */ 10.1343 + if (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360; 10.1344 + else if (lon_min > 180) lon_min -= 360 + trunc(lon_min / 360) * 360; 10.1345 + if (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360; 10.1346 + else if (lon_max > 180) lon_max -= 360 + trunc(lon_max / 360) * 360; 10.1347 + } 10.1348 + /* store rounded latitude/longitude values for round-trip safety */ 10.1349 + box->lat_min = pgl_round(lat_min); 10.1350 + box->lat_max = pgl_round(lat_max); 10.1351 + box->lon_min = pgl_round(lon_min); 10.1352 + box->lon_max = pgl_round(lon_max); 10.1353 + /* ensure that rounding does not change orientation */ 10.1354 + if (lon_min > lon_max && box->lon_min == box->lon_max) { 10.1355 + box->lon_min = -180; 10.1356 + box->lon_max = 180; 10.1357 + } 10.1358 +} 10.1359 + 10.1360 +/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */ 10.1361 +PG_FUNCTION_INFO_V1(pgl_create_ebox); 10.1362 +Datum pgl_create_ebox(PG_FUNCTION_ARGS) { 10.1363 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 10.1364 + pgl_ebox_set_boundaries( 10.1365 + box, 10.1366 + PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), 10.1367 + PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3) 10.1368 + ); 10.1369 + PG_RETURN_POINTER(box); 10.1370 +} 10.1371 + 10.1372 +/* create box ("ebox" in SQL) from two points ("epoint"s) */ 10.1373 +/* (can not be used to cover a longitude range of more than 120 degrees) */ 10.1374 +PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints); 10.1375 +Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) { 10.1376 + pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); 10.1377 + pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); 10.1378 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 10.1379 + double lat_min, lat_max, lon_min, lon_max; 10.1380 + double dlon; /* longitude range (delta longitude) */ 10.1381 + /* order latitude and longitude boundaries */ 10.1382 + if (point2->lat < point1->lat) { 10.1383 + lat_min = point2->lat; 10.1384 + lat_max = point1->lat; 10.1385 + } else { 10.1386 + lat_min = point1->lat; 10.1387 + lat_max = point2->lat; 10.1388 + } 10.1389 + if (point2->lon < point1->lon) { 10.1390 + lon_min = point2->lon; 10.1391 + lon_max = point1->lon; 10.1392 + } else { 10.1393 + lon_min = point1->lon; 10.1394 + lon_max = point2->lon; 10.1395 + } 10.1396 + /* calculate longitude range (round to avoid floating point errors) */ 10.1397 + dlon = pgl_round(lon_max - lon_min); 10.1398 + /* determine east-west direction */ 10.1399 + if (dlon >= 240) { 10.1400 + /* assume that 180th meridian is crossed and swap min/max longitude */ 10.1401 + double swap = lon_min; lon_min = lon_max; lon_max = swap; 10.1402 + } else if (dlon > 120) { 10.1403 + /* unclear orientation since delta longitude > 120 */ 10.1404 + ereport(ERROR, ( 10.1405 + errcode(ERRCODE_DATA_EXCEPTION), 10.1406 + errmsg("can not determine east/west orientation for ebox") 10.1407 + )); 10.1408 + } 10.1409 + /* use boundaries to setup box (and perform checks) */ 10.1410 + pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); 10.1411 + /* return result */ 10.1412 + PG_RETURN_POINTER(box); 10.1413 +} 10.1414 + 10.1415 +/* parse box ("ebox" in SQL) */ 10.1416 +/* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>' 10.1417 + or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */ 10.1418 +PG_FUNCTION_INFO_V1(pgl_ebox_in); 10.1419 +Datum pgl_ebox_in(PG_FUNCTION_ARGS) { 10.1420 + char *str = PG_GETARG_CSTRING(0); /* input string */ 10.1421 + char *str_lower; /* lower case version of input string */ 10.1422 + char *strptr; /* current position within string */ 10.1423 + int valid; /* number of valid chars */ 10.1424 + int done; /* specifies if latitude or longitude was read */ 10.1425 + double val; /* temporary variable */ 10.1426 + int lat_count = 0; /* count of latitude values parsed */ 10.1427 + int lon_count = 0; /* count of longitufde values parsed */ 10.1428 + double lat_min, lat_max, lon_min, lon_max; /* see pgl_box struct */ 10.1429 + pgl_box *box; /* return value (to be palloc'ed) */ 10.1430 + /* lowercase input */ 10.1431 + str_lower = psprintf("%s", str); 10.1432 + for (strptr=str_lower; *strptr; strptr++) { 10.1433 + if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; 10.1434 + } 10.1435 + /* reset reading position to start of (lowercase) string */ 10.1436 + strptr = str_lower; 10.1437 + /* check if empty box */ 10.1438 + valid = 0; 10.1439 + sscanf(strptr, " empty %n", &valid); 10.1440 + if (valid && strptr[valid] == 0) { 10.1441 + /* allocate and return empty box */ 10.1442 + box = (pgl_box *)palloc(sizeof(pgl_box)); 10.1443 + pgl_box_set_empty(box); 10.1444 + PG_RETURN_POINTER(box); 10.1445 + } 10.1446 + /* demand four blocks separated by whitespace */ 10.1447 + valid = 0; 10.1448 + sscanf(strptr, " %*s %*s %*s %*s %n", &valid); 10.1449 + /* if four blocks separated by whitespace exist, parse those blocks */ 10.1450 + if (strptr[valid] == 0) while (strptr[0]) { 10.1451 + /* parse either latitude or longitude (whichever found in input string) */ 10.1452 + done = pgl_scan(&strptr, &val, &val); 10.1453 + /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */ 10.1454 + if (done == PGL_SCAN_LAT) { 10.1455 + if (!lat_count) lat_min = val; else lat_max = val; 10.1456 + lat_count++; 10.1457 + } else if (done == PGL_SCAN_LON) { 10.1458 + if (!lon_count) lon_min = val; else lon_max = val; 10.1459 + lon_count++; 10.1460 + } else { 10.1461 + break; 10.1462 + } 10.1463 + } 10.1464 + /* require end of string, and two latitude and two longitude values */ 10.1465 + if (strptr[0] || lat_count != 2 || lon_count != 2) { 10.1466 + ereport(ERROR, ( 10.1467 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 10.1468 + errmsg("invalid input syntax for type ebox: \"%s\"", str) 10.1469 + )); 10.1470 + } 10.1471 + /* free lower case string */ 10.1472 + pfree(str_lower); 10.1473 + /* order boundaries (maximum greater than minimum) */ 10.1474 + if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; } 10.1475 + if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; } 10.1476 + /* allocate memory for result */ 10.1477 + box = (pgl_box *)palloc(sizeof(pgl_box)); 10.1478 + /* set boundaries (and perform checks) */ 10.1479 + pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max); 10.1480 + /* return result */ 10.1481 + PG_RETURN_POINTER(box); 10.1482 +} 10.1483 + 10.1484 +/* set circle to given latitude, longitude, and radius (including checks) */ 10.1485 +static void pgl_ecircle_set_latlon_radius( 10.1486 + pgl_circle *circle, double lat, double lon, double radius 10.1487 +) { 10.1488 + /* set center point (including checks) */ 10.1489 + pgl_epoint_set_latlon(&(circle->center), lat, lon); 10.1490 + /* handle non-positive radius */ 10.1491 + if (isnan(radius)) { 10.1492 + ereport(ERROR, ( 10.1493 + errcode(ERRCODE_DATA_EXCEPTION), 10.1494 + errmsg("invalid radius for ecircle") 10.1495 + )); 10.1496 + } 10.1497 + if (radius == 0) radius = 0; /* avoids -0 */ 10.1498 + else if (radius < 0) { 10.1499 + if (isfinite(radius)) { 10.1500 + ereport(NOTICE, (errmsg("negative radius converted to minus infinity"))); 10.1501 + } 10.1502 + radius = -INFINITY; 10.1503 + } 10.1504 + /* store radius (round-trip safety is ensured by pgl_print_float) */ 10.1505 + circle->radius = radius; 10.1506 +} 10.1507 + 10.1508 +/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */ 10.1509 +PG_FUNCTION_INFO_V1(pgl_create_ecircle); 10.1510 +Datum pgl_create_ecircle(PG_FUNCTION_ARGS) { 10.1511 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 10.1512 + pgl_ecircle_set_latlon_radius( 10.1513 + circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2) 10.1514 + ); 10.1515 + PG_RETURN_POINTER(circle); 10.1516 +} 10.1517 + 10.1518 +/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */ 10.1519 +PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint); 10.1520 +Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) { 10.1521 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.1522 + double radius = PG_GETARG_FLOAT8(1); 10.1523 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 10.1524 + /* set latitude, longitude, radius (and perform checks) */ 10.1525 + pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius); 10.1526 + /* return result */ 10.1527 + PG_RETURN_POINTER(circle); 10.1528 +} 10.1529 + 10.1530 +/* parse circle ("ecircle" in SQL) */ 10.1531 +/* format: '[NS]<float> [EW]<float> <float>' */ 10.1532 +PG_FUNCTION_INFO_V1(pgl_ecircle_in); 10.1533 +Datum pgl_ecircle_in(PG_FUNCTION_ARGS) { 10.1534 + char *str = PG_GETARG_CSTRING(0); /* input string */ 10.1535 + char *strptr = str; /* current position within string */ 10.1536 + double lat, lon, radius; /* parsed values as double precision flaots */ 10.1537 + int valid = 0; /* number of valid chars */ 10.1538 + int done = 0; /* stores if latitude and/or longitude was read */ 10.1539 + pgl_circle *circle; /* return value (to be palloc'ed) */ 10.1540 + /* demand three blocks separated by whitespace */ 10.1541 + sscanf(strptr, " %*s %*s %*s %n", &valid); 10.1542 + /* if three blocks separated by whitespace exist, parse those blocks */ 10.1543 + if (strptr[valid] == 0) { 10.1544 + /* parse latitude and longitude */ 10.1545 + done |= pgl_scan(&strptr, &lat, &lon); 10.1546 + done |= pgl_scan(&strptr, &lat, &lon); 10.1547 + /* parse radius (while incrementing strptr by number of bytes parsed) */ 10.1548 + valid = 0; 10.1549 + if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid; 10.1550 + } 10.1551 + /* require end of string and both latitude and longitude being parsed */ 10.1552 + if (strptr[0] || done != PGL_SCAN_LATLON) { 10.1553 + ereport(ERROR, ( 10.1554 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 10.1555 + errmsg("invalid input syntax for type ecircle: \"%s\"", str) 10.1556 + )); 10.1557 + } 10.1558 + /* allocate memory for result */ 10.1559 + circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 10.1560 + /* set latitude, longitude, radius (and perform checks) */ 10.1561 + pgl_ecircle_set_latlon_radius(circle, lat, lon, radius); 10.1562 + /* return result */ 10.1563 + PG_RETURN_POINTER(circle); 10.1564 +} 10.1565 + 10.1566 +/* parse cluster ("ecluster" in SQL) */ 10.1567 +PG_FUNCTION_INFO_V1(pgl_ecluster_in); 10.1568 +Datum pgl_ecluster_in(PG_FUNCTION_ARGS) { 10.1569 + int i; 10.1570 + char *str = PG_GETARG_CSTRING(0); /* input string */ 10.1571 + char *str_lower; /* lower case version of input string */ 10.1572 + char *strptr; /* pointer to current reading position of input */ 10.1573 + int npoints_total = 0; /* total number of points in cluster */ 10.1574 + int nentries = 0; /* total number of entries */ 10.1575 + pgl_newentry *entries; /* array of pgl_newentry to create pgl_cluster */ 10.1576 + int entries_buflen = 4; /* maximum number of elements in entries array */ 10.1577 + int valid; /* number of valid chars processed */ 10.1578 + double lat, lon; /* latitude and longitude of parsed point */ 10.1579 + int entrytype; /* current entry type */ 10.1580 + int npoints; /* number of points in current entry */ 10.1581 + pgl_point *points; /* array of pgl_point for pgl_newentry */ 10.1582 + int points_buflen; /* maximum number of elements in points array */ 10.1583 + int done; /* return value of pgl_scan function */ 10.1584 + pgl_cluster *cluster; /* created cluster */ 10.1585 + /* lowercase input */ 10.1586 + str_lower = psprintf("%s", str); 10.1587 + for (strptr=str_lower; *strptr; strptr++) { 10.1588 + if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A'; 10.1589 + } 10.1590 + /* reset reading position to start of (lowercase) string */ 10.1591 + strptr = str_lower; 10.1592 + /* allocate initial buffer for entries */ 10.1593 + entries = palloc(entries_buflen * sizeof(pgl_newentry)); 10.1594 + /* parse until end of string */ 10.1595 + while (strptr[0]) { 10.1596 + /* require previous white-space or closing parenthesis before next token */ 10.1597 + if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') { 10.1598 + goto pgl_ecluster_in_error; 10.1599 + } 10.1600 + /* ignore token "empty" */ 10.1601 + valid = 0; sscanf(strptr, " empty %n", &valid); 10.1602 + if (valid) { strptr += valid; continue; } 10.1603 + /* test for "point" token */ 10.1604 + valid = 0; sscanf(strptr, " point ( %n", &valid); 10.1605 + if (valid) { 10.1606 + strptr += valid; 10.1607 + entrytype = PGL_ENTRY_POINT; 10.1608 + goto pgl_ecluster_in_type_ok; 10.1609 + } 10.1610 + /* test for "path" token */ 10.1611 + valid = 0; sscanf(strptr, " path ( %n", &valid); 10.1612 + if (valid) { 10.1613 + strptr += valid; 10.1614 + entrytype = PGL_ENTRY_PATH; 10.1615 + goto pgl_ecluster_in_type_ok; 10.1616 + } 10.1617 + /* test for "outline" token */ 10.1618 + valid = 0; sscanf(strptr, " outline ( %n", &valid); 10.1619 + if (valid) { 10.1620 + strptr += valid; 10.1621 + entrytype = PGL_ENTRY_OUTLINE; 10.1622 + goto pgl_ecluster_in_type_ok; 10.1623 + } 10.1624 + /* test for "polygon" token */ 10.1625 + valid = 0; sscanf(strptr, " polygon ( %n", &valid); 10.1626 + if (valid) { 10.1627 + strptr += valid; 10.1628 + entrytype = PGL_ENTRY_POLYGON; 10.1629 + goto pgl_ecluster_in_type_ok; 10.1630 + } 10.1631 + /* error if no valid token found */ 10.1632 + goto pgl_ecluster_in_error; 10.1633 + pgl_ecluster_in_type_ok: 10.1634 + /* check if pgl_newentry array needs to grow */ 10.1635 + if (nentries == entries_buflen) { 10.1636 + pgl_newentry *newbuf; 10.1637 + entries_buflen *= 2; 10.1638 + newbuf = palloc(entries_buflen * sizeof(pgl_newentry)); 10.1639 + memcpy(newbuf, entries, nentries * sizeof(pgl_newentry)); 10.1640 + pfree(entries); 10.1641 + entries = newbuf; 10.1642 + } 10.1643 + /* reset number of points for current entry */ 10.1644 + npoints = 0; 10.1645 + /* allocate array for points */ 10.1646 + points_buflen = 4; 10.1647 + points = palloc(points_buflen * sizeof(pgl_point)); 10.1648 + /* parse until closing parenthesis */ 10.1649 + while (strptr[0] != ')') { 10.1650 + /* error on unexpected end of string */ 10.1651 + if (strptr[0] == 0) goto pgl_ecluster_in_error; 10.1652 + /* mark neither latitude nor longitude as read */ 10.1653 + done = PGL_SCAN_NONE; 10.1654 + /* require white-space before second, third, etc. point */ 10.1655 + if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error; 10.1656 + /* scan latitude (or longitude) */ 10.1657 + done |= pgl_scan(&strptr, &lat, &lon); 10.1658 + /* require white-space before second coordinate */ 10.1659 + if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error; 10.1660 + /* scan longitude (or latitude) */ 10.1661 + done |= pgl_scan(&strptr, &lat, &lon); 10.1662 + /* error unless both latitude and longitude were parsed */ 10.1663 + if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error; 10.1664 + /* throw error if number of points is too high */ 10.1665 + if (npoints_total == PGL_CLUSTER_MAXPOINTS) { 10.1666 + ereport(ERROR, ( 10.1667 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 10.1668 + errmsg( 10.1669 + "too many points for ecluster entry (maximum %i)", 10.1670 + PGL_CLUSTER_MAXPOINTS 10.1671 + ) 10.1672 + )); 10.1673 + } 10.1674 + /* check if pgl_point array needs to grow */ 10.1675 + if (npoints == points_buflen) { 10.1676 + pgl_point *newbuf; 10.1677 + points_buflen *= 2; 10.1678 + newbuf = palloc(points_buflen * sizeof(pgl_point)); 10.1679 + memcpy(newbuf, points, npoints * sizeof(pgl_point)); 10.1680 + pfree(points); 10.1681 + points = newbuf; 10.1682 + } 10.1683 + /* append point to pgl_point array (includes checks) */ 10.1684 + pgl_epoint_set_latlon(&(points[npoints++]), lat, lon); 10.1685 + /* increase total number of points */ 10.1686 + npoints_total++; 10.1687 + } 10.1688 + /* error if entry has no points */ 10.1689 + if (!npoints) goto pgl_ecluster_in_error; 10.1690 + /* entries with one point are automatically of type "point" */ 10.1691 + if (npoints == 1) entrytype = PGL_ENTRY_POINT; 10.1692 + /* if entries have more than one point */ 10.1693 + else { 10.1694 + /* throw error if entry type is "point" */ 10.1695 + if (entrytype == PGL_ENTRY_POINT) { 10.1696 + ereport(ERROR, ( 10.1697 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 10.1698 + errmsg("invalid input syntax for type ecluster (point entry with more than one point)") 10.1699 + )); 10.1700 + } 10.1701 + /* coerce outlines and polygons with more than 2 points to be a path */ 10.1702 + if (npoints == 2) entrytype = PGL_ENTRY_PATH; 10.1703 + } 10.1704 + /* append entry to pgl_newentry array */ 10.1705 + entries[nentries].entrytype = entrytype; 10.1706 + entries[nentries].npoints = npoints; 10.1707 + entries[nentries].points = points; 10.1708 + nentries++; 10.1709 + /* consume closing parenthesis */ 10.1710 + strptr++; 10.1711 + /* consume white-space */ 10.1712 + while (isspace(strptr[0])) strptr++; 10.1713 + } 10.1714 + /* free lower case string */ 10.1715 + pfree(str_lower); 10.1716 + /* create cluster from pgl_newentry array */ 10.1717 + cluster = pgl_new_cluster(nentries, entries); 10.1718 + /* free pgl_newentry array */ 10.1719 + for (i=0; i<nentries; i++) pfree(entries[i].points); 10.1720 + pfree(entries); 10.1721 + /* set bounding circle of cluster and check east/west orientation */ 10.1722 + if (!pgl_finalize_cluster(cluster)) { 10.1723 + ereport(ERROR, ( 10.1724 + errcode(ERRCODE_DATA_EXCEPTION), 10.1725 + errmsg("can not determine east/west orientation for ecluster"), 10.1726 + errhint("Ensure that each entry has a longitude span of less than 180 degrees.") 10.1727 + )); 10.1728 + } 10.1729 + /* return cluster */ 10.1730 + PG_RETURN_POINTER(cluster); 10.1731 + /* code to throw error */ 10.1732 + pgl_ecluster_in_error: 10.1733 + ereport(ERROR, ( 10.1734 + errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), 10.1735 + errmsg("invalid input syntax for type ecluster: \"%s\"", str) 10.1736 + )); 10.1737 +} 10.1738 + 10.1739 +/* convert point ("epoint") to string representation */ 10.1740 +PG_FUNCTION_INFO_V1(pgl_epoint_out); 10.1741 +Datum pgl_epoint_out(PG_FUNCTION_ARGS) { 10.1742 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.1743 + char latstr[PGL_NUMBUFLEN]; 10.1744 + char lonstr[PGL_NUMBUFLEN]; 10.1745 + pgl_print_lat(latstr, point->lat); 10.1746 + pgl_print_lon(lonstr, point->lon); 10.1747 + PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr)); 10.1748 +} 10.1749 + 10.1750 +/* convert box ("ebox") to string representation */ 10.1751 +PG_FUNCTION_INFO_V1(pgl_ebox_out); 10.1752 +Datum pgl_ebox_out(PG_FUNCTION_ARGS) { 10.1753 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 10.1754 + double lon_min = box->lon_min; 10.1755 + double lon_max = box->lon_max; 10.1756 + char lat_min_str[PGL_NUMBUFLEN]; 10.1757 + char lat_max_str[PGL_NUMBUFLEN]; 10.1758 + char lon_min_str[PGL_NUMBUFLEN]; 10.1759 + char lon_max_str[PGL_NUMBUFLEN]; 10.1760 + /* return string "empty" if box is set to be empty */ 10.1761 + if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty"); 10.1762 + /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */ 10.1763 + /* (required since pgl_box_in orders the longitude boundaries) */ 10.1764 + if (lon_min > lon_max) { 10.1765 + if (lon_min + lon_max >= 0) lon_min -= 360; 10.1766 + else lon_max += 360; 10.1767 + } 10.1768 + /* format and return result */ 10.1769 + pgl_print_lat(lat_min_str, box->lat_min); 10.1770 + pgl_print_lat(lat_max_str, box->lat_max); 10.1771 + pgl_print_lon(lon_min_str, lon_min); 10.1772 + pgl_print_lon(lon_max_str, lon_max); 10.1773 + PG_RETURN_CSTRING(psprintf( 10.1774 + "%s %s %s %s", 10.1775 + lat_min_str, lon_min_str, lat_max_str, lon_max_str 10.1776 + )); 10.1777 +} 10.1778 + 10.1779 +/* convert circle ("ecircle") to string representation */ 10.1780 +PG_FUNCTION_INFO_V1(pgl_ecircle_out); 10.1781 +Datum pgl_ecircle_out(PG_FUNCTION_ARGS) { 10.1782 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 10.1783 + char latstr[PGL_NUMBUFLEN]; 10.1784 + char lonstr[PGL_NUMBUFLEN]; 10.1785 + char radstr[PGL_NUMBUFLEN]; 10.1786 + pgl_print_lat(latstr, circle->center.lat); 10.1787 + pgl_print_lon(lonstr, circle->center.lon); 10.1788 + pgl_print_float(radstr, circle->radius); 10.1789 + PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr)); 10.1790 +} 10.1791 + 10.1792 +/* convert cluster ("ecluster") to string representation */ 10.1793 +PG_FUNCTION_INFO_V1(pgl_ecluster_out); 10.1794 +Datum pgl_ecluster_out(PG_FUNCTION_ARGS) { 10.1795 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 10.1796 + char latstr[PGL_NUMBUFLEN]; /* string buffer for latitude */ 10.1797 + char lonstr[PGL_NUMBUFLEN]; /* string buffer for longitude */ 10.1798 + char ***strings; /* array of array of strings */ 10.1799 + char *string; /* string of current token */ 10.1800 + char *res, *resptr; /* result and pointer to current write position */ 10.1801 + size_t reslen = 1; /* length of result (init with 1 for terminator) */ 10.1802 + int npoints; /* number of points of current entry */ 10.1803 + int i, j; /* i: entry, j: point in entry */ 10.1804 + /* handle empty clusters */ 10.1805 + if (cluster->nentries == 0) { 10.1806 + /* free detoasted cluster (if copy) */ 10.1807 + PG_FREE_IF_COPY(cluster, 0); 10.1808 + /* return static result */ 10.1809 + PG_RETURN_CSTRING("empty"); 10.1810 + } 10.1811 + /* allocate array of array of strings */ 10.1812 + strings = palloc(cluster->nentries * sizeof(char **)); 10.1813 + /* iterate over all entries in cluster */ 10.1814 + for (i=0; i<cluster->nentries; i++) { 10.1815 + /* get number of points in entry */ 10.1816 + npoints = cluster->entries[i].npoints; 10.1817 + /* allocate array of strings (one string for each point plus two extra) */ 10.1818 + strings[i] = palloc((2 + npoints) * sizeof(char *)); 10.1819 + /* determine opening string */ 10.1820 + switch (cluster->entries[i].entrytype) { 10.1821 + case PGL_ENTRY_POINT: string = (i==0)?"point (" :" point ("; break; 10.1822 + case PGL_ENTRY_PATH: string = (i==0)?"path (" :" path ("; break; 10.1823 + case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break; 10.1824 + case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break; 10.1825 + default: string = (i==0)?"unknown" :" unknown"; 10.1826 + } 10.1827 + /* use opening string as first string in array */ 10.1828 + strings[i][0] = string; 10.1829 + /* update result length (for allocating result string later) */ 10.1830 + reslen += strlen(string); 10.1831 + /* iterate over all points */ 10.1832 + for (j=0; j<npoints; j++) { 10.1833 + /* create string representation of point */ 10.1834 + pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat); 10.1835 + pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon); 10.1836 + string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr); 10.1837 + /* copy string pointer to string array */ 10.1838 + strings[i][j+1] = string; 10.1839 + /* update result length (for allocating result string later) */ 10.1840 + reslen += strlen(string); 10.1841 + } 10.1842 + /* use closing parenthesis as last string in array */ 10.1843 + strings[i][npoints+1] = ")"; 10.1844 + /* update result length (for allocating result string later) */ 10.1845 + reslen++; 10.1846 + } 10.1847 + /* allocate result string */ 10.1848 + res = palloc(reslen); 10.1849 + /* set write pointer to begin of result string */ 10.1850 + resptr = res; 10.1851 + /* copy strings into result string */ 10.1852 + for (i=0; i<cluster->nentries; i++) { 10.1853 + npoints = cluster->entries[i].npoints; 10.1854 + for (j=0; j<npoints+2; j++) { 10.1855 + string = strings[i][j]; 10.1856 + strcpy(resptr, string); 10.1857 + resptr += strlen(string); 10.1858 + /* free strings allocated by psprintf */ 10.1859 + if (j != 0 && j != npoints+1) pfree(string); 10.1860 + } 10.1861 + /* free array of strings */ 10.1862 + pfree(strings[i]); 10.1863 + } 10.1864 + /* free array of array of strings */ 10.1865 + pfree(strings); 10.1866 + /* free detoasted cluster (if copy) */ 10.1867 + PG_FREE_IF_COPY(cluster, 0); 10.1868 + /* return result */ 10.1869 + PG_RETURN_CSTRING(res); 10.1870 +} 10.1871 + 10.1872 +/* binary input function for point ("epoint") */ 10.1873 +PG_FUNCTION_INFO_V1(pgl_epoint_recv); 10.1874 +Datum pgl_epoint_recv(PG_FUNCTION_ARGS) { 10.1875 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 10.1876 + pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point)); 10.1877 + point->lat = pq_getmsgfloat8(buf); 10.1878 + point->lon = pq_getmsgfloat8(buf); 10.1879 + PG_RETURN_POINTER(point); 10.1880 +} 10.1881 + 10.1882 +/* binary input function for box ("ebox") */ 10.1883 +PG_FUNCTION_INFO_V1(pgl_ebox_recv); 10.1884 +Datum pgl_ebox_recv(PG_FUNCTION_ARGS) { 10.1885 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 10.1886 + pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box)); 10.1887 + box->lat_min = pq_getmsgfloat8(buf); 10.1888 + box->lat_max = pq_getmsgfloat8(buf); 10.1889 + box->lon_min = pq_getmsgfloat8(buf); 10.1890 + box->lon_max = pq_getmsgfloat8(buf); 10.1891 + PG_RETURN_POINTER(box); 10.1892 +} 10.1893 + 10.1894 +/* binary input function for circle ("ecircle") */ 10.1895 +PG_FUNCTION_INFO_V1(pgl_ecircle_recv); 10.1896 +Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) { 10.1897 + StringInfo buf = (StringInfo)PG_GETARG_POINTER(0); 10.1898 + pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle)); 10.1899 + circle->center.lat = pq_getmsgfloat8(buf); 10.1900 + circle->center.lon = pq_getmsgfloat8(buf); 10.1901 + circle->radius = pq_getmsgfloat8(buf); 10.1902 + PG_RETURN_POINTER(circle); 10.1903 +} 10.1904 + 10.1905 +/* TODO: binary receive function for cluster */ 10.1906 + 10.1907 +/* binary output function for point ("epoint") */ 10.1908 +PG_FUNCTION_INFO_V1(pgl_epoint_send); 10.1909 +Datum pgl_epoint_send(PG_FUNCTION_ARGS) { 10.1910 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.1911 + StringInfoData buf; 10.1912 + pq_begintypsend(&buf); 10.1913 + pq_sendfloat8(&buf, point->lat); 10.1914 + pq_sendfloat8(&buf, point->lon); 10.1915 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 10.1916 +} 10.1917 + 10.1918 +/* binary output function for box ("ebox") */ 10.1919 +PG_FUNCTION_INFO_V1(pgl_ebox_send); 10.1920 +Datum pgl_ebox_send(PG_FUNCTION_ARGS) { 10.1921 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 10.1922 + StringInfoData buf; 10.1923 + pq_begintypsend(&buf); 10.1924 + pq_sendfloat8(&buf, box->lat_min); 10.1925 + pq_sendfloat8(&buf, box->lat_max); 10.1926 + pq_sendfloat8(&buf, box->lon_min); 10.1927 + pq_sendfloat8(&buf, box->lon_max); 10.1928 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 10.1929 +} 10.1930 + 10.1931 +/* binary output function for circle ("ecircle") */ 10.1932 +PG_FUNCTION_INFO_V1(pgl_ecircle_send); 10.1933 +Datum pgl_ecircle_send(PG_FUNCTION_ARGS) { 10.1934 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 10.1935 + StringInfoData buf; 10.1936 + pq_begintypsend(&buf); 10.1937 + pq_sendfloat8(&buf, circle->center.lat); 10.1938 + pq_sendfloat8(&buf, circle->center.lon); 10.1939 + pq_sendfloat8(&buf, circle->radius); 10.1940 + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); 10.1941 +} 10.1942 + 10.1943 +/* TODO: binary send functions for cluster */ 10.1944 + 10.1945 +/* cast point ("epoint") to box ("ebox") */ 10.1946 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox); 10.1947 +Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) { 10.1948 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.1949 + pgl_box *box = palloc(sizeof(pgl_box)); 10.1950 + box->lat_min = point->lat; 10.1951 + box->lat_max = point->lat; 10.1952 + box->lon_min = point->lon; 10.1953 + box->lon_max = point->lon; 10.1954 + PG_RETURN_POINTER(box); 10.1955 +} 10.1956 + 10.1957 +/* cast point ("epoint") to circle ("ecircle") */ 10.1958 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle); 10.1959 +Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) { 10.1960 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.1961 + pgl_circle *circle = palloc(sizeof(pgl_box)); 10.1962 + circle->center = *point; 10.1963 + circle->radius = 0; 10.1964 + PG_RETURN_POINTER(circle); 10.1965 +} 10.1966 + 10.1967 +/* cast point ("epoint") to cluster ("ecluster") */ 10.1968 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster); 10.1969 +Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) { 10.1970 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.1971 + pgl_newentry entry; 10.1972 + entry.entrytype = PGL_ENTRY_POINT; 10.1973 + entry.npoints = 1; 10.1974 + entry.points = point; 10.1975 + PG_RETURN_POINTER(pgl_new_cluster(1, &entry)); 10.1976 +} 10.1977 + 10.1978 +/* cast box ("ebox") to cluster ("ecluster") */ 10.1979 +#define pgl_ebox_to_ecluster_macro(i, a, b) \ 10.1980 + entries[i].entrytype = PGL_ENTRY_POLYGON; \ 10.1981 + entries[i].npoints = 4; \ 10.1982 + entries[i].points = points[i]; \ 10.1983 + points[i][0].lat = box->lat_min; \ 10.1984 + points[i][0].lon = (a); \ 10.1985 + points[i][1].lat = box->lat_min; \ 10.1986 + points[i][1].lon = (b); \ 10.1987 + points[i][2].lat = box->lat_max; \ 10.1988 + points[i][2].lon = (b); \ 10.1989 + points[i][3].lat = box->lat_max; \ 10.1990 + points[i][3].lon = (a); 10.1991 +PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster); 10.1992 +Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) { 10.1993 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0); 10.1994 + double lon, dlon; 10.1995 + int nentries; 10.1996 + pgl_newentry entries[3]; 10.1997 + pgl_point points[3][4]; 10.1998 + if (box->lat_min > box->lat_max) { 10.1999 + nentries = 0; 10.2000 + } else if (box->lon_min > box->lon_max) { 10.2001 + if (box->lon_min < 0) { 10.2002 + lon = pgl_round((box->lon_min + 180) / 2.0); 10.2003 + nentries = 3; 10.2004 + pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); 10.2005 + pgl_ebox_to_ecluster_macro(1, lon, 180); 10.2006 + pgl_ebox_to_ecluster_macro(2, -180, box->lon_max); 10.2007 + } else if (box->lon_max > 0) { 10.2008 + lon = pgl_round((box->lon_max - 180) / 2.0); 10.2009 + nentries = 3; 10.2010 + pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); 10.2011 + pgl_ebox_to_ecluster_macro(1, -180, lon); 10.2012 + pgl_ebox_to_ecluster_macro(2, lon, box->lon_max); 10.2013 + } else { 10.2014 + nentries = 2; 10.2015 + pgl_ebox_to_ecluster_macro(0, box->lon_min, 180); 10.2016 + pgl_ebox_to_ecluster_macro(1, -180, box->lon_max); 10.2017 + } 10.2018 + } else { 10.2019 + dlon = pgl_round(box->lon_max - box->lon_min); 10.2020 + if (dlon < 180) { 10.2021 + nentries = 1; 10.2022 + pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max); 10.2023 + } else { 10.2024 + lon = pgl_round((box->lon_min + box->lon_max) / 2.0); 10.2025 + if ( 10.2026 + pgl_round(lon - box->lon_min) < 180 && 10.2027 + pgl_round(box->lon_max - lon) < 180 10.2028 + ) { 10.2029 + nentries = 2; 10.2030 + pgl_ebox_to_ecluster_macro(0, box->lon_min, lon); 10.2031 + pgl_ebox_to_ecluster_macro(1, lon, box->lon_max); 10.2032 + } else { 10.2033 + nentries = 3; 10.2034 + pgl_ebox_to_ecluster_macro(0, box->lon_min, -60); 10.2035 + pgl_ebox_to_ecluster_macro(1, -60, 60); 10.2036 + pgl_ebox_to_ecluster_macro(2, 60, box->lon_max); 10.2037 + } 10.2038 + } 10.2039 + } 10.2040 + PG_RETURN_POINTER(pgl_new_cluster(nentries, entries)); 10.2041 +} 10.2042 + 10.2043 +/* extract latitude from point ("epoint") */ 10.2044 +PG_FUNCTION_INFO_V1(pgl_epoint_lat); 10.2045 +Datum pgl_epoint_lat(PG_FUNCTION_ARGS) { 10.2046 + PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat); 10.2047 +} 10.2048 + 10.2049 +/* extract longitude from point ("epoint") */ 10.2050 +PG_FUNCTION_INFO_V1(pgl_epoint_lon); 10.2051 +Datum pgl_epoint_lon(PG_FUNCTION_ARGS) { 10.2052 + PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon); 10.2053 +} 10.2054 + 10.2055 +/* extract minimum latitude from box ("ebox") */ 10.2056 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_min); 10.2057 +Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) { 10.2058 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min); 10.2059 +} 10.2060 + 10.2061 +/* extract maximum latitude from box ("ebox") */ 10.2062 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_max); 10.2063 +Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) { 10.2064 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max); 10.2065 +} 10.2066 + 10.2067 +/* extract minimum longitude from box ("ebox") */ 10.2068 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_min); 10.2069 +Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) { 10.2070 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min); 10.2071 +} 10.2072 + 10.2073 +/* extract maximum longitude from box ("ebox") */ 10.2074 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_max); 10.2075 +Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) { 10.2076 + PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max); 10.2077 +} 10.2078 + 10.2079 +/* extract center point from circle ("ecircle") */ 10.2080 +PG_FUNCTION_INFO_V1(pgl_ecircle_center); 10.2081 +Datum pgl_ecircle_center(PG_FUNCTION_ARGS) { 10.2082 + PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center)); 10.2083 +} 10.2084 + 10.2085 +/* extract radius from circle ("ecircle") */ 10.2086 +PG_FUNCTION_INFO_V1(pgl_ecircle_radius); 10.2087 +Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) { 10.2088 + PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius); 10.2089 +} 10.2090 + 10.2091 +/* check if point is inside box (overlap operator "&&") in SQL */ 10.2092 +PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap); 10.2093 +Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) { 10.2094 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.2095 + pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1); 10.2096 + PG_RETURN_BOOL(pgl_point_in_box(point, box)); 10.2097 +} 10.2098 + 10.2099 +/* check if point is inside circle (overlap operator "&&") in SQL */ 10.2100 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap); 10.2101 +Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) { 10.2102 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.2103 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 10.2104 + PG_RETURN_BOOL( 10.2105 + pgl_distance( 10.2106 + point->lat, point->lon, 10.2107 + circle->center.lat, circle->center.lon 10.2108 + ) <= circle->radius 10.2109 + ); 10.2110 +} 10.2111 + 10.2112 +/* check if point is inside cluster (overlap operator "&&") in SQL */ 10.2113 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap); 10.2114 +Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) { 10.2115 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.2116 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 10.2117 + bool retval = pgl_point_in_cluster(point, cluster); 10.2118 + PG_FREE_IF_COPY(cluster, 1); 10.2119 + PG_RETURN_BOOL(retval); 10.2120 +} 10.2121 + 10.2122 +/* check if two boxes overlap (overlap operator "&&") in SQL */ 10.2123 +PG_FUNCTION_INFO_V1(pgl_ebox_overlap); 10.2124 +Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) { 10.2125 + pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0); 10.2126 + pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1); 10.2127 + PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2)); 10.2128 +} 10.2129 + 10.2130 +/* check if two circles overlap (overlap operator "&&") in SQL */ 10.2131 +PG_FUNCTION_INFO_V1(pgl_ecircle_overlap); 10.2132 +Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) { 10.2133 + pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); 10.2134 + pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); 10.2135 + PG_RETURN_BOOL( 10.2136 + pgl_distance( 10.2137 + circle1->center.lat, circle1->center.lon, 10.2138 + circle2->center.lat, circle2->center.lon 10.2139 + ) <= circle1->radius + circle2->radius 10.2140 + ); 10.2141 +} 10.2142 + 10.2143 +/* check if circle and cluster overlap (overlap operator "&&") in SQL */ 10.2144 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap); 10.2145 +Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) { 10.2146 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 10.2147 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 10.2148 + bool retval = ( 10.2149 + pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius 10.2150 + ); 10.2151 + PG_FREE_IF_COPY(cluster, 1); 10.2152 + PG_RETURN_BOOL(retval); 10.2153 +} 10.2154 + 10.2155 +/* calculate distance between two points ("<->" operator) in SQL */ 10.2156 +PG_FUNCTION_INFO_V1(pgl_epoint_distance); 10.2157 +Datum pgl_epoint_distance(PG_FUNCTION_ARGS) { 10.2158 + pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0); 10.2159 + pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1); 10.2160 + PG_RETURN_FLOAT8(pgl_distance( 10.2161 + point1->lat, point1->lon, point2->lat, point2->lon 10.2162 + )); 10.2163 +} 10.2164 + 10.2165 +/* calculate point to circle distance ("<->" operator) in SQL */ 10.2166 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance); 10.2167 +Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) { 10.2168 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.2169 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1); 10.2170 + double distance = pgl_distance( 10.2171 + point->lat, point->lon, circle->center.lat, circle->center.lon 10.2172 + ) - circle->radius; 10.2173 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 10.2174 +} 10.2175 + 10.2176 +/* calculate point to cluster distance ("<->" operator) in SQL */ 10.2177 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance); 10.2178 +Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) { 10.2179 + pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0); 10.2180 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 10.2181 + double distance = pgl_point_cluster_distance(point, cluster); 10.2182 + PG_FREE_IF_COPY(cluster, 1); 10.2183 + PG_RETURN_FLOAT8(distance); 10.2184 +} 10.2185 + 10.2186 +/* calculate distance between two circles ("<->" operator) in SQL */ 10.2187 +PG_FUNCTION_INFO_V1(pgl_ecircle_distance); 10.2188 +Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) { 10.2189 + pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0); 10.2190 + pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1); 10.2191 + double distance = pgl_distance( 10.2192 + circle1->center.lat, circle1->center.lon, 10.2193 + circle2->center.lat, circle2->center.lon 10.2194 + ) - (circle1->radius + circle2->radius); 10.2195 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 10.2196 +} 10.2197 + 10.2198 +/* calculate circle to cluster distance ("<->" operator) in SQL */ 10.2199 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance); 10.2200 +Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) { 10.2201 + pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0); 10.2202 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 10.2203 + double distance = ( 10.2204 + pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius 10.2205 + ); 10.2206 + PG_FREE_IF_COPY(cluster, 1); 10.2207 + PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance); 10.2208 +} 10.2209 + 10.2210 + 10.2211 +/*-----------------------------------------------------------* 10.2212 + * B-tree comparison operators and index support functions * 10.2213 + *-----------------------------------------------------------*/ 10.2214 + 10.2215 +/* macro for a B-tree operator (without detoasting) */ 10.2216 +#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \ 10.2217 + PG_FUNCTION_INFO_V1(func); \ 10.2218 + Datum func(PG_FUNCTION_ARGS) { \ 10.2219 + type *a = (type *)PG_GETARG_POINTER(0); \ 10.2220 + type *b = (type *)PG_GETARG_POINTER(1); \ 10.2221 + PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \ 10.2222 + } 10.2223 + 10.2224 +/* macro for a B-tree comparison function (without detoasting) */ 10.2225 +#define PGL_BTREE_CMP(func, type, cmpfunc) \ 10.2226 + PG_FUNCTION_INFO_V1(func); \ 10.2227 + Datum func(PG_FUNCTION_ARGS) { \ 10.2228 + type *a = (type *)PG_GETARG_POINTER(0); \ 10.2229 + type *b = (type *)PG_GETARG_POINTER(1); \ 10.2230 + PG_RETURN_INT32(cmpfunc(a, b)); \ 10.2231 + } 10.2232 + 10.2233 +/* macro for a B-tree operator (with detoasting) */ 10.2234 +#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \ 10.2235 + PG_FUNCTION_INFO_V1(func); \ 10.2236 + Datum func(PG_FUNCTION_ARGS) { \ 10.2237 + bool res; \ 10.2238 + type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ 10.2239 + type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ 10.2240 + res = cmpfunc(a, b) oper 0; \ 10.2241 + PG_FREE_IF_COPY(a, 0); \ 10.2242 + PG_FREE_IF_COPY(b, 1); \ 10.2243 + PG_RETURN_BOOL(res); \ 10.2244 + } 10.2245 + 10.2246 +/* macro for a B-tree comparison function (with detoasting) */ 10.2247 +#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \ 10.2248 + PG_FUNCTION_INFO_V1(func); \ 10.2249 + Datum func(PG_FUNCTION_ARGS) { \ 10.2250 + int32_t res; \ 10.2251 + type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \ 10.2252 + type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \ 10.2253 + res = cmpfunc(a, b); \ 10.2254 + PG_FREE_IF_COPY(a, 0); \ 10.2255 + PG_FREE_IF_COPY(b, 1); \ 10.2256 + PG_RETURN_INT32(res); \ 10.2257 + } 10.2258 + 10.2259 +/* B-tree operators and comparison function for point */ 10.2260 +PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <) 10.2261 +PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=) 10.2262 +PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==) 10.2263 +PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=) 10.2264 +PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=) 10.2265 +PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >) 10.2266 +PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp) 10.2267 + 10.2268 +/* B-tree operators and comparison function for box */ 10.2269 +PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <) 10.2270 +PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=) 10.2271 +PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==) 10.2272 +PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=) 10.2273 +PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=) 10.2274 +PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >) 10.2275 +PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp) 10.2276 + 10.2277 +/* B-tree operators and comparison function for circle */ 10.2278 +PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <) 10.2279 +PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=) 10.2280 +PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==) 10.2281 +PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=) 10.2282 +PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=) 10.2283 +PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >) 10.2284 +PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp) 10.2285 + 10.2286 + 10.2287 +/*--------------------------------* 10.2288 + * GiST index support functions * 10.2289 + *--------------------------------*/ 10.2290 + 10.2291 +/* GiST "consistent" support function */ 10.2292 +PG_FUNCTION_INFO_V1(pgl_gist_consistent); 10.2293 +Datum pgl_gist_consistent(PG_FUNCTION_ARGS) { 10.2294 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 10.2295 + pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); 10.2296 + StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); 10.2297 + bool *recheck = (bool *)PG_GETARG_POINTER(4); 10.2298 + /* demand recheck because index and query methods are lossy */ 10.2299 + *recheck = true; 10.2300 + /* strategy number 11: equality of two points */ 10.2301 + if (strategy == 11) { 10.2302 + /* query datum is another point */ 10.2303 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 10.2304 + /* convert other point to key */ 10.2305 + pgl_pointkey querykey; 10.2306 + pgl_point_to_key(query, querykey); 10.2307 + /* return true if both keys overlap */ 10.2308 + PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); 10.2309 + } 10.2310 + /* strategy number 13: equality of two circles */ 10.2311 + if (strategy == 13) { 10.2312 + /* query datum is another circle */ 10.2313 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 10.2314 + /* convert other circle to key */ 10.2315 + pgl_areakey querykey; 10.2316 + pgl_circle_to_key(query, querykey); 10.2317 + /* return true if both keys overlap */ 10.2318 + PG_RETURN_BOOL(pgl_keys_overlap(key, querykey)); 10.2319 + } 10.2320 + /* for all remaining strategies, keys on empty objects produce no match */ 10.2321 + /* (check necessary because query radius may be infinite) */ 10.2322 + if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false); 10.2323 + /* strategy number 21: overlapping with point */ 10.2324 + if (strategy == 21) { 10.2325 + /* query datum is a point */ 10.2326 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 10.2327 + /* return true if estimated distance (allowed to be smaller than real 10.2328 + distance) between index key and point is zero */ 10.2329 + PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0); 10.2330 + } 10.2331 + /* strategy number 22: (point) overlapping with box */ 10.2332 + if (strategy == 22) { 10.2333 + /* query datum is a box */ 10.2334 + pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1); 10.2335 + /* determine bounding box of indexed key */ 10.2336 + pgl_box keybox; 10.2337 + pgl_key_to_box(key, &keybox); 10.2338 + /* return true if query box overlaps with bounding box of indexed key */ 10.2339 + PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox)); 10.2340 + } 10.2341 + /* strategy number 23: overlapping with circle */ 10.2342 + if (strategy == 23) { 10.2343 + /* query datum is a circle */ 10.2344 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 10.2345 + /* return true if estimated distance (allowed to be smaller than real 10.2346 + distance) between index key and circle center is smaller than radius */ 10.2347 + PG_RETURN_BOOL( 10.2348 + pgl_estimate_key_distance(key, &(query->center)) <= query->radius 10.2349 + ); 10.2350 + } 10.2351 + /* strategy number 24: overlapping with cluster */ 10.2352 + if (strategy == 24) { 10.2353 + bool retval; /* return value */ 10.2354 + /* query datum is a cluster */ 10.2355 + pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 10.2356 + /* return true if estimated distance (allowed to be smaller than real 10.2357 + distance) between index key and circle center is smaller than radius */ 10.2358 + retval = ( 10.2359 + pgl_estimate_key_distance(key, &(query->bounding.center)) <= 10.2360 + query->bounding.radius 10.2361 + ); 10.2362 + PG_FREE_IF_COPY(query, 1); /* free detoasted cluster (if copy) */ 10.2363 + PG_RETURN_BOOL(retval); 10.2364 + } 10.2365 + /* throw error for any unknown strategy number */ 10.2366 + elog(ERROR, "unrecognized strategy number: %d", strategy); 10.2367 +} 10.2368 + 10.2369 +/* GiST "union" support function */ 10.2370 +PG_FUNCTION_INFO_V1(pgl_gist_union); 10.2371 +Datum pgl_gist_union(PG_FUNCTION_ARGS) { 10.2372 + GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); 10.2373 + pgl_keyptr out; /* return value (to be palloc'ed) */ 10.2374 + int i; 10.2375 + /* determine key size */ 10.2376 + size_t keysize = PGL_KEY_IS_AREAKEY( 10.2377 + (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key) 10.2378 + ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey); 10.2379 + /* begin with first key as result */ 10.2380 + out = palloc(keysize); 10.2381 + memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize); 10.2382 + /* unite current result with second, third, etc. key */ 10.2383 + for (i=1; i<entryvec->n; i++) { 10.2384 + pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)); 10.2385 + } 10.2386 + /* return result */ 10.2387 + PG_RETURN_POINTER(out); 10.2388 +} 10.2389 + 10.2390 +/* GiST "compress" support function for indicis on points */ 10.2391 +PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint); 10.2392 +Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) { 10.2393 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 10.2394 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 10.2395 + /* only transform new leaves */ 10.2396 + if (entry->leafkey) { 10.2397 + /* get point to be transformed */ 10.2398 + pgl_point *point = (pgl_point *)DatumGetPointer(entry->key); 10.2399 + /* allocate memory for key */ 10.2400 + pgl_keyptr key = palloc(sizeof(pgl_pointkey)); 10.2401 + /* transform point to key */ 10.2402 + pgl_point_to_key(point, key); 10.2403 + /* create new GISTENTRY structure as return value */ 10.2404 + retval = palloc(sizeof(GISTENTRY)); 10.2405 + gistentryinit( 10.2406 + *retval, PointerGetDatum(key), 10.2407 + entry->rel, entry->page, entry->offset, FALSE 10.2408 + ); 10.2409 + } else { 10.2410 + /* inner nodes have already been transformed */ 10.2411 + retval = entry; 10.2412 + } 10.2413 + /* return pointer to old or new GISTENTRY structure */ 10.2414 + PG_RETURN_POINTER(retval); 10.2415 +} 10.2416 + 10.2417 +/* GiST "compress" support function for indicis on circles */ 10.2418 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle); 10.2419 +Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) { 10.2420 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 10.2421 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 10.2422 + /* only transform new leaves */ 10.2423 + if (entry->leafkey) { 10.2424 + /* get circle to be transformed */ 10.2425 + pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key); 10.2426 + /* allocate memory for key */ 10.2427 + pgl_keyptr key = palloc(sizeof(pgl_areakey)); 10.2428 + /* transform circle to key */ 10.2429 + pgl_circle_to_key(circle, key); 10.2430 + /* create new GISTENTRY structure as return value */ 10.2431 + retval = palloc(sizeof(GISTENTRY)); 10.2432 + gistentryinit( 10.2433 + *retval, PointerGetDatum(key), 10.2434 + entry->rel, entry->page, entry->offset, FALSE 10.2435 + ); 10.2436 + } else { 10.2437 + /* inner nodes have already been transformed */ 10.2438 + retval = entry; 10.2439 + } 10.2440 + /* return pointer to old or new GISTENTRY structure */ 10.2441 + PG_RETURN_POINTER(retval); 10.2442 +} 10.2443 + 10.2444 +/* GiST "compress" support function for indices on clusters */ 10.2445 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster); 10.2446 +Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) { 10.2447 + GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); 10.2448 + GISTENTRY *retval; /* return value (to be palloc'ed unless set to entry) */ 10.2449 + /* only transform new leaves */ 10.2450 + if (entry->leafkey) { 10.2451 + /* get cluster to be transformed (detoasting necessary!) */ 10.2452 + pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key); 10.2453 + /* allocate memory for key */ 10.2454 + pgl_keyptr key = palloc(sizeof(pgl_areakey)); 10.2455 + /* transform cluster to key */ 10.2456 + pgl_circle_to_key(&(cluster->bounding), key); 10.2457 + /* create new GISTENTRY structure as return value */ 10.2458 + retval = palloc(sizeof(GISTENTRY)); 10.2459 + gistentryinit( 10.2460 + *retval, PointerGetDatum(key), 10.2461 + entry->rel, entry->page, entry->offset, FALSE 10.2462 + ); 10.2463 + /* free detoasted datum */ 10.2464 + if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster); 10.2465 + } else { 10.2466 + /* inner nodes have already been transformed */ 10.2467 + retval = entry; 10.2468 + } 10.2469 + /* return pointer to old or new GISTENTRY structure */ 10.2470 + PG_RETURN_POINTER(retval); 10.2471 +} 10.2472 + 10.2473 +/* GiST "decompress" support function for indices */ 10.2474 +PG_FUNCTION_INFO_V1(pgl_gist_decompress); 10.2475 +Datum pgl_gist_decompress(PG_FUNCTION_ARGS) { 10.2476 + /* return passed pointer without transformation */ 10.2477 + PG_RETURN_POINTER(PG_GETARG_POINTER(0)); 10.2478 +} 10.2479 + 10.2480 +/* GiST "penalty" support function */ 10.2481 +PG_FUNCTION_INFO_V1(pgl_gist_penalty); 10.2482 +Datum pgl_gist_penalty(PG_FUNCTION_ARGS) { 10.2483 + GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0); 10.2484 + GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1); 10.2485 + float *penalty = (float *)PG_GETARG_POINTER(2); 10.2486 + /* get original key and key to insert */ 10.2487 + pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key); 10.2488 + pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key); 10.2489 + /* copy original key */ 10.2490 + union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key; 10.2491 + if (PGL_KEY_IS_AREAKEY(orig)) { 10.2492 + memcpy(union_key.areakey, orig, sizeof(union_key.areakey)); 10.2493 + } else { 10.2494 + memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey)); 10.2495 + } 10.2496 + /* calculate union of both keys */ 10.2497 + pgl_unite_keys((pgl_keyptr)&union_key, new); 10.2498 + /* penalty equal to reduction of key length (logarithm of added area) */ 10.2499 + /* (return value by setting referenced value and returning pointer) */ 10.2500 + *penalty = ( 10.2501 + PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key) 10.2502 + ); 10.2503 + PG_RETURN_POINTER(penalty); 10.2504 +} 10.2505 + 10.2506 +/* GiST "picksplit" support function */ 10.2507 +PG_FUNCTION_INFO_V1(pgl_gist_picksplit); 10.2508 +Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) { 10.2509 + GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0); 10.2510 + GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1); 10.2511 + OffsetNumber i; /* between FirstOffsetNumber and entryvec->n (inclusive) */ 10.2512 + union { 10.2513 + pgl_pointkey pointkey; 10.2514 + pgl_areakey areakey; 10.2515 + } union_all; /* union of all keys (to be calculated from scratch) 10.2516 + (later cut in half) */ 10.2517 + int is_areakey = PGL_KEY_IS_AREAKEY( 10.2518 + (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key) 10.2519 + ); 10.2520 + int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey); 10.2521 + pgl_keyptr unionL = palloc(keysize); /* union of keys that go left */ 10.2522 + pgl_keyptr unionR = palloc(keysize); /* union of keys that go right */ 10.2523 + pgl_keyptr key; /* current key to be processed */ 10.2524 + /* allocate memory for array of left and right keys, set counts to zero */ 10.2525 + v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); 10.2526 + v->spl_nleft = 0; 10.2527 + v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber)); 10.2528 + v->spl_nright = 0; 10.2529 + /* calculate union of all keys from scratch */ 10.2530 + memcpy( 10.2531 + (pgl_keyptr)&union_all, 10.2532 + (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key), 10.2533 + keysize 10.2534 + ); 10.2535 + for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) { 10.2536 + pgl_unite_keys( 10.2537 + (pgl_keyptr)&union_all, 10.2538 + (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key) 10.2539 + ); 10.2540 + } 10.2541 + /* check if trivial split is necessary due to exhausted key length */ 10.2542 + /* (Note: keys for empty objects must have node depth set to maximum) */ 10.2543 + if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == ( 10.2544 + is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH 10.2545 + )) { 10.2546 + /* half of all keys go left */ 10.2547 + for ( 10.2548 + i=FirstOffsetNumber; 10.2549 + i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; 10.2550 + i=OffsetNumberNext(i) 10.2551 + ) { 10.2552 + /* pointer to current key */ 10.2553 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 10.2554 + /* update unionL */ 10.2555 + /* check if key is first key that goes left */ 10.2556 + if (!v->spl_nleft) { 10.2557 + /* first key that goes left is just copied to unionL */ 10.2558 + memcpy(unionL, key, keysize); 10.2559 + } else { 10.2560 + /* unite current value and next key */ 10.2561 + pgl_unite_keys(unionL, key); 10.2562 + } 10.2563 + /* append offset number to list of keys that go left */ 10.2564 + v->spl_left[v->spl_nleft++] = i; 10.2565 + } 10.2566 + /* other half goes right */ 10.2567 + for ( 10.2568 + i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2; 10.2569 + i<entryvec->n; 10.2570 + i=OffsetNumberNext(i) 10.2571 + ) { 10.2572 + /* pointer to current key */ 10.2573 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 10.2574 + /* update unionR */ 10.2575 + /* check if key is first key that goes right */ 10.2576 + if (!v->spl_nright) { 10.2577 + /* first key that goes right is just copied to unionR */ 10.2578 + memcpy(unionR, key, keysize); 10.2579 + } else { 10.2580 + /* unite current value and next key */ 10.2581 + pgl_unite_keys(unionR, key); 10.2582 + } 10.2583 + /* append offset number to list of keys that go right */ 10.2584 + v->spl_right[v->spl_nright++] = i; 10.2585 + } 10.2586 + } 10.2587 + /* otherwise, a non-trivial split is possible */ 10.2588 + else { 10.2589 + /* cut covered area in half */ 10.2590 + /* (union_all then refers to area of keys that go left) */ 10.2591 + /* check if union of all keys covers empty and non-empty objects */ 10.2592 + if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) { 10.2593 + /* if yes, split into empty and non-empty objects */ 10.2594 + pgl_key_set_empty((pgl_keyptr)&union_all); 10.2595 + } else { 10.2596 + /* otherwise split by next bit */ 10.2597 + ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++; 10.2598 + /* NOTE: type bit conserved */ 10.2599 + } 10.2600 + /* determine for each key if it goes left or right */ 10.2601 + for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) { 10.2602 + /* pointer to current key */ 10.2603 + key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key); 10.2604 + /* keys within one half of the area go left */ 10.2605 + if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) { 10.2606 + /* update unionL */ 10.2607 + /* check if key is first key that goes left */ 10.2608 + if (!v->spl_nleft) { 10.2609 + /* first key that goes left is just copied to unionL */ 10.2610 + memcpy(unionL, key, keysize); 10.2611 + } else { 10.2612 + /* unite current value of unionL and processed key */ 10.2613 + pgl_unite_keys(unionL, key); 10.2614 + } 10.2615 + /* append offset number to list of keys that go left */ 10.2616 + v->spl_left[v->spl_nleft++] = i; 10.2617 + } 10.2618 + /* the other keys go right */ 10.2619 + else { 10.2620 + /* update unionR */ 10.2621 + /* check if key is first key that goes right */ 10.2622 + if (!v->spl_nright) { 10.2623 + /* first key that goes right is just copied to unionR */ 10.2624 + memcpy(unionR, key, keysize); 10.2625 + } else { 10.2626 + /* unite current value of unionR and processed key */ 10.2627 + pgl_unite_keys(unionR, key); 10.2628 + } 10.2629 + /* append offset number to list of keys that go right */ 10.2630 + v->spl_right[v->spl_nright++] = i; 10.2631 + } 10.2632 + } 10.2633 + } 10.2634 + /* store unions in return value */ 10.2635 + v->spl_ldatum = PointerGetDatum(unionL); 10.2636 + v->spl_rdatum = PointerGetDatum(unionR); 10.2637 + /* return all results */ 10.2638 + PG_RETURN_POINTER(v); 10.2639 +} 10.2640 + 10.2641 +/* GiST "same"/"equal" support function */ 10.2642 +PG_FUNCTION_INFO_V1(pgl_gist_same); 10.2643 +Datum pgl_gist_same(PG_FUNCTION_ARGS) { 10.2644 + pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0); 10.2645 + pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1); 10.2646 + bool *boolptr = (bool *)PG_GETARG_POINTER(2); 10.2647 + /* two keys are equal if they are binary equal */ 10.2648 + /* (return result by setting referenced boolean and returning pointer) */ 10.2649 + *boolptr = !memcmp( 10.2650 + key1, 10.2651 + key2, 10.2652 + PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey) 10.2653 + ); 10.2654 + PG_RETURN_POINTER(boolptr); 10.2655 +} 10.2656 + 10.2657 +/* GiST "distance" support function */ 10.2658 +PG_FUNCTION_INFO_V1(pgl_gist_distance); 10.2659 +Datum pgl_gist_distance(PG_FUNCTION_ARGS) { 10.2660 + GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0); 10.2661 + pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key); 10.2662 + StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2); 10.2663 + bool *recheck = (bool *)PG_GETARG_POINTER(4); 10.2664 + double distance; /* return value */ 10.2665 + /* demand recheck because distance is just an estimation */ 10.2666 + /* (real distance may be bigger) */ 10.2667 + *recheck = true; 10.2668 + /* strategy number 31: distance to point */ 10.2669 + if (strategy == 31) { 10.2670 + /* query datum is a point */ 10.2671 + pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1); 10.2672 + /* use pgl_estimate_pointkey_distance() function to compute result */ 10.2673 + distance = pgl_estimate_key_distance(key, query); 10.2674 + /* avoid infinity (reserved!) */ 10.2675 + if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 10.2676 + /* return result */ 10.2677 + PG_RETURN_FLOAT8(distance); 10.2678 + } 10.2679 + /* strategy number 33: distance to circle */ 10.2680 + if (strategy == 33) { 10.2681 + /* query datum is a circle */ 10.2682 + pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1); 10.2683 + /* estimate distance to circle center and substract circle radius */ 10.2684 + distance = ( 10.2685 + pgl_estimate_key_distance(key, &(query->center)) - query->radius 10.2686 + ); 10.2687 + /* convert non-positive values to zero and avoid infinity (reserved!) */ 10.2688 + if (distance <= 0) distance = 0; 10.2689 + else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 10.2690 + /* return result */ 10.2691 + PG_RETURN_FLOAT8(distance); 10.2692 + } 10.2693 + /* strategy number 34: distance to cluster */ 10.2694 + if (strategy == 34) { 10.2695 + /* query datum is a cluster */ 10.2696 + pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 10.2697 + /* estimate distance to bounding center and substract bounding radius */ 10.2698 + distance = ( 10.2699 + pgl_estimate_key_distance(key, &(query->bounding.center)) - 10.2700 + query->bounding.radius 10.2701 + ); 10.2702 + /* convert non-positive values to zero and avoid infinity (reserved!) */ 10.2703 + if (distance <= 0) distance = 0; 10.2704 + else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE; 10.2705 + /* free detoasted cluster (if copy) */ 10.2706 + PG_FREE_IF_COPY(query, 1); 10.2707 + /* return result */ 10.2708 + PG_RETURN_FLOAT8(distance); 10.2709 + } 10.2710 + /* throw error for any unknown strategy number */ 10.2711 + elog(ERROR, "unrecognized strategy number: %d", strategy); 10.2712 +} 10.2713 +
11.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 11.2 +++ b/latlon.control Sun Aug 21 17:43:48 2016 +0200 11.3 @@ -0,0 +1,4 @@ 11.4 +# geographical extension 11.5 +comment = 'geospatial support' 11.6 +default_version = '0.1' 11.7 +relocatable = true
12.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 12.2 +++ b/make-doc.sh Sun Aug 21 17:43:48 2016 +0200 12.3 @@ -0,0 +1,8 @@ 12.4 +#!/bin/sh 12.5 +# 12.6 +# This command can be used to update the README.html file after changing the 12.7 +# README.mkd file. 12.8 + 12.9 +echo "<html><head><title>"`grep '[^ \t\r\n][^ \t\r\n]*' README.mkd | head -n 1`"</title></head><body>" > README.html 12.10 +markdown2 README.mkd >> README.html 12.11 +echo "</body></html>" >> README.html