The upcoming jOOQ 3.16 will finally offer support for the various RDBMS GIS extensions via issue #982. This is great news per se, and will be covered in a future blog post, when the integration is ready. This post here is about something else.
Adding support for such a feature is a great source of procrastination. You see, when you develop a backend library, all you ever do is implement algorithms, API, and other abstract stuff. But GIS is different. With GIS, a SQL developer like me suddenly has access to a “UI”, and having access to a “UI” woke the inner programmer child in me.
I was pleasantly surprised that DBeaver, my preferred SQL editor, has out of the box support for GIS’s WKT format. E.g. run a query like this:
SELECT st_geomfromtext('polygon((0 0, 2 3, 0 6, -2 3, 0 0))');
The output being
Let’s play around with GIS
So, naturally, what other thing to do than play around with it? The most obvious next thing to generate would be your favourite logo, which happens to be very easy to map to polygons. Let’s look at a few GIS features step by step. And for this blog post, I’ll be using PostGIS, which you can get very easily from docker hub:
WITH
-- These "sprites" are just reusable 2D polygons, which I declared
-- to avoid having to deal with the lengthy WKT format
sprites AS (
SELECT
-- We project the original 1x1 square, and a scaled 4x4 version
s AS square,
st_scale(s, 4, 4) AS square4
-- Here, we're creating a square called "s"
FROM (VALUES
(st_polygonfromtext('polygon ((0 0, 1 0, 1 1, 0 1, 0 0))'))
) t (s)
)
-- st_union combines 2 polygons into a multipolygon
SELECT st_union(square, st_translate(square4, 2, 0))
FROM sprites
The output of the above query is
MULTIPOLYGON (((0 0, 1 0, 1 1, 0 1, 0 0)), ((2 0, 6 0, 6 4, 2 4, 2 0)))
Or, using DBeaver’s visualisation utility:
The st_union isn’t much different from any other set union. Note that I’ve translated the larger square 2 units to the right, so they don’t overlap. Otherwise, the union would have just been the bigger square.
A polygon describes a set of points in the 2D number plane. Creating a union of those two sets of points is natural. We can also create a difference of the two squares, instead:
WITH
sprites AS (
SELECT
s AS square,
st_scale(s, 4, 4) AS square4
FROM (VALUES
(st_polygonfromtext('polygon ((0 0, 1 0, 1 1, 0 1, 0 0))'))
) t (s)
)
SELECT st_difference(square4, square)
FROM sprites
Resulting in:
POLYGON ((0 4, 4 4, 4 0, 1 0, 1 1, 0 1, 0 4))
Or, visually:
With these simple tools, we can now create all 4 letters of the jOOQ logo. As a reminder, the tools were:
st_polygonfromtext: Create a polygon from a WKT text representationst_scale: Scale any geometry to a new sizest_translate: Translate any geometry to a new positionst_union: Combine two geometries (or more, if used as an aggregate function)st_difference: Remove one geometry from another
The GIS API is vast, both in PostGIS as well as in the ISO/IEC 13249-3:2016 standard version, but these simple tools shall suffice for now. Let’s look at this query:
WITH
-- Creating the two squares to work with
sprites AS (
SELECT
s AS square,
st_scale(s, 4, 4) AS square4
FROM (VALUES
(st_polygonfromtext('polygon ((0 0, 1 0, 1 1, 0 1, 0 0))'))
) t (s)
),
-- All the 4 letters j, o1, o2, q of the jOOQ logo
letters AS (
SELECT
-- The letter j
st_difference(
st_difference(
st_difference(
square4,
st_translate(st_scale(square, 2, 3), 1, 1)
),
st_translate(st_scale(square, 1, 2), 0, 2)
),
st_translate(st_scale(square, 1, 0.5), 3, 2.5)
) AS j,
-- The letter o
st_difference(square4, st_translate(square, 1, 1)) AS o1,
-- The letter o
st_difference(square4, st_translate(square, 2, 2)) AS o2,
-- The letter q
st_union(
st_difference(
square4,
st_translate(st_scale(square, 2, 2), 1, 1)
),
st_translate(st_scale(square, 1, 2.5), 1.5, -1)
) as q
FROM sprites
)
SELECT
-- Combine all the 4 letters
st_union(v)
FROM
letters,
-- Arrange the 4 letters next to each other
LATERAL (VALUES
(st_translate(j, 0, 5)),
(st_translate(o1, 5, 5)),
(o2),
(st_translate(q, 5, 0))
) t (v);
This produces:
Neat, huh?
Next step: The Mandelbrot Set
A natural next step for any procrastinator is to generate the Mandelbrot set. To prepare ourselves for what’s behind the Mandelbrot, have a look at this neat video (more procrastination):
There are different ways to calculate it with SQL, here’s one that I came up with:
WITH RECURSIVE
-- These are the dimensions that you can play around with
dims (r1, r2, i1, i2, s, it, p) AS (
VALUES (
-- The dimensions of the real axis
-2::float, 1::float,
-- The dimensions of the imaginary axis
-1.5::float, 1.5::float,
-- The step size on each axis, per pixel or sprite
0.01::float,
-- The maximum number of iterations
100,
-- "Infinity", i.e. when to stop
256.0::float
)
),
-- The square again, as before
sprites (s) AS (VALUES
(st_polygonfromtext('polygon ((0 0, 0 1, 1 1, 1 0, 0 0))'))
),
-- The number plane, as ints
n1 (r, i) AS (
SELECT r, i
FROM
dims,
generate_series((r1 / s)::int, (r2 / s)::int) r,
generate_series((i1 / s)::int, (i2 / s)::int) i
),
-- The number plane as scaled floats
n2 (r, i) AS (
SELECT r::float * s::float, i::float * s::float
FROM dims, n1
),
-- The recursive calculation of the Mandelbrot formula
-- zn = (zn-1)^2 + c
l (cr, ci, zr, zi, g, it, p) AS (
SELECT r, i, 0::float, 0::float, 0, it, p FROM n2, dims
UNION ALL
SELECT cr, ci, zr*zr - zi*zi + cr, 2*zr*zi + ci, g + 1, it, p
FROM l
-- The recursions stops when we reach the maximum
WHERE g < it
-- Or, when we reach "infinity"
AND zr*zr + zi*zi < p
),
-- Find the last calculated value per point in the
-- complex number plane c (cr, ci), discard the others
x (cr, ci, zr, zi, g) AS (
SELECT DISTINCT ON (cr, ci) cr, ci, zr, zi, g
FROM l
ORDER BY cr, ci, g DESC
)
-- Turn every calculated point into a square
SELECT
st_union(
st_translate(sprites.s, round(cr / dims.s), round(ci / dims.s))
)
FROM x, sprites, dims
-- Retain only the points *not* belonging to the Mandelbrot set
-- You can also inverse the equation to retain the points that
-- belong to the set
WHERE zr*zr + zi*zi > p;
Note, I’m using an artificial “infinity” here, because:
- That speeds things up with little noticeable difference at this zoom scale
- I couldn’t figure out how to make PostgreSQL overflow float operations to float infinities, like Java or CockroachDB or other IEEE 754 float implementations do. Any help appreciated, see this Stack Overflow question.
The output is the well known shape
You can play around with this and use different values for “dims” to zoom in, e.g.
dims (r1, r2, i1, i2, s, it, p) as (values (
(-0.925-0.032)::float, (-0.925+0.032)::float,
(0.266-0.032)::float, (0.266+0.032)::float,
0.00005::float, 100, 256.0::float)
)
This will generate some really neat zooms, all with SQL:
Granted, it’s not the most efficient way to calculate these things, but that wasn’t the point here, was it?
Game of Life
I can be nerd-sniped:
And of course, I had to accept that challenge too! Now, the Game of Life is a simple “game” where we have the x,y natural number plane (e.g. “pixels”), and with each coordinate, we associate whether that “cell” is dead or alive. Then, we establish the following set of rules:
- Any live cell with two or three live neighbours survives.
- Any dead cell with three live neighbours becomes a live cell.
- All other live cells die in the next generation. Similarly, all other dead cells stay dead.
No, it’s hard to “animate” things in SQL using spatial extensions, so as a workaround, I’ll just display iterations of 100×100 pixel tiles of the Game of Life next to each other. The first iteration is just the starting point, e.g. a random number of “alive” cells:
WITH RECURSIVE
-- Again generate a "sprite"a from a square polygon
sprites (s) AS (VALUES (
st_polygonfromtext('polygon ((0 0, 0 1, 1 1, 1 0, 0 0))')
)),
-- Generate the x, y number plane and associate a boolean value
-- with each coordinate, generated randomly.
-- 10% of all cells are alive
m (x, y, b) AS (
SELECT x, y,
CASE WHEN random() > 0.9 THEN 1 ELSE 0 END
FROM
generate_series(1, 100) x,
generate_series(1, 100) y
)
SELECT st_union(st_translate(sprites.s, x::float, y::float))
FROM m, sprites
-- Display only "alive" cells
WHERE b = 1
This will produce something like:
Now, all we have to do is iterate the game formula. Now, recursive SQL has quite a few limitations. E.g. I couldn’t get PostgreSQL to aggregate recursive data, nor self-join the recursive table to find the nearest neighbors of any cell. But with window functions, this is definitely possible. So here goes:
WITH RECURSIVE
sprites (s) AS (VALUES (
st_polygonfromtext('polygon ((0 0, 0 1, 1 1, 1 0, 0 0))')
)),
m (x, y, b) AS (
SELECT x, y,
CASE WHEN random() > 0.9 THEN 1 ELSE 0 END
FROM
generate_series(1, 100) x,
generate_series(1, 100) y
),
-- This is the recursion of the Game of Life
e (x, y, b, g) AS (
-- The recursion starts with the above initial data
SELECT x, y, b, 1 FROM m
UNION ALL
SELECT
e.x, e.y,
CASE
-- If this cell is alive, and 2 or 3
-- neighbors are alive, too, then this
-- cell stays alive
WHEN e.b = 1 AND
sum(e.b) OVER w1
+ sum(e.b) OVER w2
+ sum(e.b) OVER w3 IN (2, 3)
THEN 1
-- If this cell is dead, and 3 neighbors
-- are alive, then this cell becomes alive
WHEN e.b = 0 AND
sum(e.b) OVER w1
+ sum(e.b) OVER w2
+ sum(e.b) OVER w3 = 3
THEN 1
-- Otherwise, this cell dies or stays dead
ELSE 0
END,
e.g + 1
FROM e
-- The recursion stops after 100 iterations
WHERE e.g < 100
WINDOW
-- We order the data set by x, y. In SQL, there isn't
-- a notion of a 2-dimensional number plane. We only
-- have 1 dimension.
o AS (ORDER BY x, y),
-- w1 is all the neighbors of the previous row in the x, y
-- plane, which is all the SQL rows that are 101-99 rows
-- before the current row.
w1 AS (o ROWS BETWEEN 101 PRECEDING AND 99 PRECEDING),
-- w2 is all the neighbors of the current row in the x, y
-- number plane, excluding the current cell
w2 AS (o ROWS BETWEEN 1 PRECEDING AND 1 FOLLOWING
EXCLUDE CURRENT ROW),
-- w3 is all the neighbors of the next row
w3 AS (o ROWS BETWEEN 99 FOLLOWING AND 101 FOLLOWING)
)
-- Finally, combine all the iterations
SELECT st_union(st_translate(
sprites.s,
(x + (g - 1) % 10 * 120)::float,
(y - (g - 1) / 10 * 120)::float
))
FROM e, sprites
WHERE b = 1
;
Think of the window specifications as follows:
----+-------------+-------------+-------------+-------------+---- ... | 105: (1, 5) | 106: (1, 6) | 107: (1, 7) | 108: (1, 8) | ... ... | 205: (2, 5) | 206: (2, 6) | 207: (2, 7) | 208: (2, 8) | ... ... | 305: (3, 5) | 306: (3, 6) | 307: (3, 7) | 308: (3, 8) | ... ... | 405: (4, 5) | 406: (4, 6) | 407: (4, 7) | 408: (4, 8) | ...
So, in a 100×100 grid, x=3,y=7 can be encoded as 307, and its neighbors are
- w1: 206, 207, 208, i.e. between 101 preceding and 99 preceding of 307
- w2: 306, 308, i.e. between 1 preceding and 1 following of 307
- w3: 406, 407, 408, i.e. between 99 following and 101 following of 307
The output looks like this:
Or, zooming in on iterations 1-4:
Or 21-24:
That’s really cool, huh!
Glider Gun
The nerd-snipe was to animate the glider gun pattern, which just means we’ll have to replace the random generation of the first iteration by something constant.
WITH RECURSIVE
sprites (s) AS (VALUES (
st_polygonfromtext('polygon ((0 0, 0 1, 1 1, 1 0, 0 0))')
)),
m (x, y, b) AS (
SELECT x, y,
-- Initial data here
CASE WHEN (x, y) IN (
(2, 6), (2, 7), (3, 6), (3, 7), (12, 6), (12, 7), (12, 8),
(13, 5), (13, 9), (14, 4), (14, 10), (15, 4), (15, 10),
(16, 7), (17, 5), (17, 9), (18, 6), (18, 7), (18, 8),
(19, 7), (22, 4), (22, 5), (22, 6), (23, 4), (23, 5),
(23, 6), (24, 3), (24, 7), (26, 2), (26, 3), (26, 7),
(26, 8), (36, 4), (36, 5), (37, 4), (37, 5)
) THEN 1 ELSE 0 END
FROM
generate_series(1, 100) x,
generate_series(1, 100) y
),
e (x, y, b, g) AS (
SELECT x, y, b, 1 FROM m
UNION ALL
SELECT
e.x, e.y,
CASE
WHEN e.b = 1 AND
sum(e.b) OVER w1
+ sum(e.b) OVER w2
+ sum(e.b) OVER w3 IN (2, 3)
THEN 1
WHEN e.b = 0 AND
sum(e.b) OVER w1
+ sum(e.b) OVER w2
+ sum(e.b) OVER w3 = 3
THEN 1
ELSE 0
END,
e.g + 1
FROM e
WHERE e.g < 100
WINDOW
o AS (ORDER BY x, y),
w1 AS (o ROWS BETWEEN 101 PRECEDING AND 99 PRECEDING),
w2 AS (o ROWS BETWEEN 1 PRECEDING AND 1 FOLLOWING
EXCLUDE CURRENT ROW),
w3 AS (o ROWS BETWEEN 99 FOLLOWING AND 101 FOLLOWING)
)
SELECT st_union(st_translate(
sprites.s,
(x + (g - 1) % 10 * 120)::float,
(y - (g - 1) / 10 * 120)::float
))
FROM e, sprites
WHERE b = 1
;
And, as you can see, the gun works:
It started with this:
And eventually produces the well known gliders:
from Java, SQL and jOOQ. https://ift.tt/322vhw2
via IFTTT
No comments:
Post a Comment