2022 / Day 30: Remix
I was really planning to rework 2022 / Day 19: Globe to x-axis meters in logarithmic scale but it looked really boring. So instead today I’m going to do a 1-degree graticule (re: 2022 / Day 13: A five minute map) for a flat earth map (re: 2022 / Day 20: Favourite).
The SQL itself is much simpler because it essentially only needs to draw a set
of concentric rings around the Null Island spanning out 180 degrees with a 1
degree step (as parallels). This can be done straight on the location, not
needing to st_translate everything
from the pole (because there’s no need to measure any real distances). Meaning
an incremental st_buffer coupled
with st_boundary will suffice. This
is done in the CTE paral.
For the meridian lines from the Null Island the CTE called merid constructs
a line from the origin to every point measured 1 degree apart along a circle
that covers the whole World - essentially a
doing a st_lineinterpolatepoint
every 1/360 fraction along the
st_boundary of a
st_buffer set at 180 degrees from
the point of origin st_point(0,0,4326)::geometry. It’s important to do this
using the geometry type, not geography because otherwise the Earth would be
assumed to be round. Also, as the first node a buffer will be at 90 degrees
clockwise from north, the query includes a
st_rotate
for radians(90) so the first point will be at
st_azimuth 0.
The transformation from epsg:4326 to
epsg:3857 is done using QGIS. Because PostGIS will
complain with a transform: tolerance condition error (-20). QGIS will survive
with linestrings going out of the epsg:3857 bounds (
too far North/South) it seems as long as the first and last point of the linestring
are within the bounds. That’s the reason why the merid lines are chopped to
the correct length with doing a
st_intersection with the World’s
bounds from the CTE bounds.
with
bounds as (
select
st_envelope(
st_collect(
st_point(-180,89,4326),
st_point(180,-89,4326)
)
) as geom
),
paral as (
select
s, s::varchar as lbl,
case
when mod(s,5)=0 then 1
else 0
end as deg,
case
when mod(s,10)=0 then 1
else 0
end as as_lbl,
st_boundary(
st_buffer(
st_point(0, 0, 4326),
s,
'quad_segs=80'
)
) geom
from
generate_series(0, 180, 1)s
),
merid as (
select
s,
case
/* don't show labels for these*/
when s = 0 then null
when s = 90 then null
/* otherwise go with pi fractions*/
when s = 10 then 'π/18'
when s = 20 then 'π/9'
when s = 30 then 'π/6'
when s = 40 then '2π/9'
when s = 50 then '5π/18'
when s = 60 then 'π/3'
when s = 70 then '7π/18'
when s = 80 then '4π/9'
when s = 90 then 'π/2'
when s = 100 then '5π/9'
when s = 110 then '11π/18'
when s = 120 then '2π/3'
when s = 130 then '13π/18'
when s = 140 then '7π/9'
when s = 150 then '5π/6'
when s = 160 then '8π/9'
when s = 170 then '17π/18'
when s = 180 then 'π'
when s = 190 then '-17π/18'
when s = 200 then '-8π/9'
when s = 210 then '-5π/6'
when s = 220 then '-7π/9'
when s = 230 then '-13π/18'
when s = 240 then '-2π/3'
when s = 250 then '-11π/18'
when s = 260 then '-5π/9'
when s = 270 then '-π/2'
when s = 280 then '-4π/9'
when s = 290 then '-7π/18'
when s = 300 then '-π/3'
when s = 310 then '-5π/18'
when s = 320 then '-2π/9)'
when s = 330 then '-π/6)'
when s = 340 then '-π/9'
when s = 350 then '-π/18'
when s = 360 then '-π'
end as lbl,
case
when mod(s,5)=0 then 1
else 0
end as deg,
case
when mod(s,10)=0 then 1
else 0
end as as_lbl,
st_intersection(b.geom, d.geom) as geom
from (
select
s,
st_makeline(
st_point(0,0,4326),
st_lineinterpolatepoint(
st_boundary(
st_rotate(
st_buffer(
st_point(0,0,4326),
180
),
radians(90)
)
),
s/360.0
)
) as geom
from
generate_series(0,360,1)s
) d, bounds b
)
select
row_number() over()::int as oid,
s, lbl, deg, as_lbl, cl, geom
from (
select s, lbl, deg, as_lbl, geom, 'merid' as cl from merid
union all
select s, lbl, deg, as_lbl, geom, 'paral' as cl from paral
) d
;