2022 / Day 21: Kontur World Population
Todays query will group the world population at the Equator and then draw a joyline based on the population values derived from the Kontur World Population dataset. There might be other reasons but mainly I’m doing this so I could pack all the data into the query below (as is the tradition with all my maps for this years 30DayMapChallenge).
Anyway. For data preprocessing I imported the
Kontur World Population
dataset using ogr2ogr to my local
PostgreSQL database with:
$ ogr2ogr -f \
PostgreSQL \
-nlt GEOMETRY \
"PG:host=localhost user=postgres dbname=postgres" \
/home/data/kontur_population.gpkg
And then in the database to preapre the data for todays query:
select
st_asencodedpolyline(
st_makeline(
array_agg(
st_translate(
st_force2d(pts.geom),
0,
st_z(pts.geom)/1000000.0
) order by pts.path
)
) as geom
from (
select
1 as oid,
st_segmentize(
st_makeline(
array_agg(
st_pointz(
q_lon,
0,
population,
4326
) order by q_lon
)
),
1
) as geom
from (
select
sum(coalesce(population, 0)) as population, q_lon
from (
select
st_x(
st_snaptogrid(
st_transform(
st_centroid(geom),
4326
),
1
)
) as q_lon, population
from
population
) f
group by q_lon
order by q_lon
) f
) d
join lateral
st_dumppoints(d.geom) pts on true
group by
d.oid;
Which aggregates the h3 hex population values into 1 degree bins using
st_snaptogrid. Then aggregates
these into a 3D line along the Equator with the sum of population as
z-coordinate. st_segmentize
the resulting line at 1 degree, dumps the points in a lateral join and
then st_translate every nodes
y-coordinate by a value of population / 1000000.0. And finally outputs the
linestring in polyline format using
st_asencodedpolyline. This
is the string we’ll use as an input for todays querys data.
So we’ll start the CTE with decoding the geometry with st_linefromencodedpolyline. We’ll plot the data on an ellipse centered at Null Island and the semi-minor axis reaching the poles and semi-major axis along the Equator as far as possible. This will give the result an oblique view centered at around [0,45]. As yesterday we’ll create a full set of 360 pieces of them rotating every one by 1 degree with st_rotate. We’ll need the boundary only (and hence st_boundary). And another thing here - we’ll reverse the resulting linestrings with st_reverse because the outer shell of a polygon runs clockwise, but we’ll want it counter-clockwise :)
From the population data in pop.geom we’ll extract all vertices at all
timestamps we have for ellipses, find the corresponding location on the ellipse
boundary with st_lineinterpolatepoint
and st_translate the nodes y-coordinate
the required amount base on the st_y value
from the population-polyline
With to give a bit of context let’s add the longitudinal locations of the 20 biggest cities we were looking at a few days ago in 2022 / Day 19: Globe as well on the Equator line. The list is compiled according to Wikipedias [List of largest cities](https://en.wikipedia.org/wiki/List_of_largest_cities with their locations from NaturalEarth.
The process here is similar only that we will not have to shift any y-coordinate
values - only st_lineinterpolatepoint
the correct amount.
with
pop as (
select
st_linefromencodedpolyline(
'{pA~fsia@be@_ibEhL_ibEbR_ibEuv@_ibE{|M_ibEbxK_ibEha@_ibEw`d@_ibE`yZ_ibEnzJ_ibEtM_ibEnA_ibEeG_ibEgA_ibEbB_ibExE_ibEqC_ibEsA_ibErB_ibEa`G_ibEioA_ibEawxD_ibE~nwD_ibEe~]_ibE`jQ_ibEhoW_ibEsB_ibEc_D_ibEewF_ibEuulA_ibEbjz@_ibEdgY_ibEjqB_ibEkV_ibEhj@_ibE|M_ibEzH_ibEgC_ibEgI_ibEoP_ibErB_ibEjX_ibEY_ibE_L_ibEeoJ_ibEdqH_ibE~e@_ibEunA_ibEptB_ibEibB_ibEiqA_ibEn{B_ibEyzB_ibEpqA_ibEaj\_ibEqjdC_ibE{dzn@_ibEgqeX_ibEdu}o@_ibEz_jH_ibEmsjD_ibEqqjz@_ibEhwzS_ibEdhds@_ibEymaJ_ibEfi`A_ibElofK_ibEewel@_ibE|`y_@_ibE|cnF_ibEbubC_ibEcmL_ibEcf{J_ibEosjE_ibEevoU_ibE|sj_@_ibEu|_j@_ibE`c_G_ibEwvqJ_ibEyjut@_ibEizjhA_ibEhchmA_ibEc_yD_ibErk~z@_ibEcrhA_ibE`moV_ibEe_}Y_ibEs]_ibEqisO_ibEuqhL_ibExqeW_ibE}`pu@_ibEjrbp@_ibEongL_ibEb`hZ_ibEw}yr@_ibEjo|]_ibEeaaV_ibErgdF_ibEewfvA_ibEfgtw@_ibEjwzp@_ibE}~j|B_ibElm|s@_ibE|g}Z_ibEejroB_ibEnwdmB_ibEh_qC_ibEyosZ_ibE|hlfA_ibEt`dY_ibEnlY_ibEkbmD_ibEreaP_ibE}pfD_ibEa`mB_ibEjtwL_ibExgeM_ibEgh~S_ibEjs`J_ibEityc@_ibEm|rN_ibEzkd|@_ibEstxQ_ibEt~uA_ibE|ngP_ibEu|mB_ibEiw}K_ibEk_o[_ibEhpcV_ibEaq``A_ibEnonJ_ibEo}kdA_ibEj~|zA_ibErf}d@_ibEmetr@_ibEehhU_ibE|f~dA_ibExydE_ibEkilJ_ibEe_}O_ibElgeK_ibEx_eS_ibEs|{E_ibEyxnR_ibEzxmv@_ibEm@_ibEm@_ibEu`@_ibEpc@_ibEizB_ibEel@_ibEijE_ibE{wJ_ibEggM_ibEirY_ibErxf@_ibEaeU_ibEzmi@_ibEno@_ibEdl@_ibE_|i@_ibEa~nt@_ibEo|vC_ibEnqtc@_ibEemgJ_ibEcwmG_ibElkfH_ibEouxE_ibEpu~E_ibEce~cA_ibEq~ab@_ibErm`b@_ibEiy`i@_ibEnvgf@_ibEg`tdA_ibEz{ud@_ibE{fmxA_ibErmuD_ibEulh]_ibE|c}jB_ibEa{s~@_ibEyumbA_ibEvawfA_ibEaxtkA_ibExz~D_ibEqf{zC_ibElnnk@_ibEfql`@_ibE~kugA_ibE`|_~@_ibEj~kN_ibEu~x{A_ibE~_q_B_ibE||I_ibEb{e{@_ibEjyxX_ibEmilK_ibE_dsh@_ibElt}h@_ibEwuwT_ibEbvdS_ibElsl@_ibEeona@_ibEj~cY_ibEodyH_ibEezax@_ibEsv}aA_ibEqjnuB_ibEwczY_ibEwfuuG_ibEzwp}J_ibEu~n~@_ibEfjvwA_ibEorhzA_ibExhja@_ibEsihq@_ibEmmvH_ibE~`|d@_ibEr{{|B_ibEvqc|B_ibEdahJ_ibEcigu@_ibEugraC_ibEhlkoB_ibE`alq@_ibEgp|h@_ibEtrv@_ibEfg`|@_ibEwhE_ibE}pek@_ibE`rndA_ibE`~fY_ibEjkuJ_ibEgz_T_ibEmbyE_ibE`zh[_ibE}jlK_ibE|{vO_ibEejqP_ibEjtO_ibErs{N_ibExouN_ibE_|_J_ibEokeF_ibE}erM_ibE{`ruB_ibEvfwU_ibEmefrA_ibEzucfA_ibEopm~B_ibEc~qrA_ibEupcyK_ibEbohgC_ibEj{ekA_ibEo`m{A_ibEe|urF_ibEzffqA_ibEladyE_ibE{vwW_ibEztopE_ibEhwjK_ibE`{tJ_ibEtyzj@_ibEyhwaB_ibE`rxE_ibEdnwcB_ibEc}rxD_ibErw~fD_ibE_unr@_ibEr~veB_ibExkexA_ibEr_fbC_ibEjerH_ibEk_zi@_ibE}{j_@_ibEtzyyA_ibEgxxS_ibEsfbv@_ibEsyv_B_ibEsu~[_ibEtkz^_ibE_x`U_ibEkkbuD_ibEuoy]_ibEubkhB_ibEsymxA_ibE~`pkG_ibEpcuC_ibE}b|\_ibEts}Y_ibE}nwcB_ibEiealG_ibEdbbS_ibEbimtC_ibEugwmB_ibE`fqn@_ibEr{dxD_ibE_k}G_ibEepijB_ibEwimuF_ibEdon_O_ibEneo_B_ibEljzP_ibEgwxQ_ibEjomZ_ibEcjxlC_ibEjkm{F_ibEyvyN_ibEll}R_ibExcmD_ibEnvcg@_ibE}wH_ibEzu}B_ibE}`hd@_ibEo{jN_ibEsbeA_ibE|jun@_ibEirueA_ibEwcdoB_ibEzik`D_ibEdwhl@_ibE|fR_ibEgtlB_ibEiq|W_ibEfqa^_ibEmoj@_ibEn~uG_ibEsrk@_ibE~eY_ibEuvha@_ibEjr~\_ibEuwyM_ibE|myU_ibE{bM_ibEbg^_ibEwcD_ibEm`D_ibE}kP_ibEg|J_ibEpfO_ibEhkZ_ibEllE_ibEd_@_ibEceC_ibEs{d@_ibElnP_ibEunZ_ibEz}[_ibEzaC_ibE{wO_ibEfpL_ibEingC_ibEhsyA_ibEgzrL_ibErk~J_ibE|tf@_ibEinL_ibEh}q@_ibEdri@_ibE'
) geom
),
cities as (
select
name, cast(lon as numeric) lon,
cast(lat as numeric),
round(cast(lat as numeric),1) as rounded_lat
from (
values (
'Tokyo', 139.749462, 35.686963
),(
'Delhi', 77.199980, 28.600023
),(
'Shanghai', 121.434559, 31.218398
),(
'São Paulo', -46.626966, -23.556734
),(
'Mexico City', -99.132934, 19.444388
),(
'Cairo', 31.248022, 30.051906
),(
'Mumbai', 72.855043, 19.018936
),(
'Beijing', 116.386340, 39.930838
),(
'Dhaka', 90.406634, 23.725006
),(
'Osaka', 135.458199, 34.751981
),(
'New York', -73.981963, 40.751925
),(
'Karachi', 67.01, 24.86
),(
'Buenos Aires', -58.399477, -34.600556
),(
'Chongqing', 106.5504, 29.5637
),(
'Istanbul', 29.008056, 41.106942
),(
'Kolkata',88.37, 22.5675
),(
'Manila',120.980271, 14.606105
),(
'Lagos',3.389585,6.445208
),(
'Rio do Janeiro',-43.226967, -22.923077
),(
'Tianjin', 117.2054, 39.1336
)
) x (name, lon, lat)
),
ellipsis as (
select
s,
st_reverse(
st_boundary(
st_scale(
st_rotate(
st_buffer(
st_point(0,0, 4326),
180.0
),
radians(s),
st_point(0,0,4326)
),
1, 0.5
)
)
) as geom
from
generate_series(0,360, 1) s
),
pop_lines as (
select
s,
st_setsrid(
st_makeline(
array_agg(
st_translate(
st_force2d(pts),
0,
st_y(pops.geom)
) order by p
)
),
4326) as geom
from
ellipsis d
join lateral
generate_series(0,360,1) p on true
join lateral
st_lineinterpolatepoint(
d.geom,
p::numeric/360.0
) pts on true
left join (
select
p.path[1]-1 as path, p.geom
from
pop
join lateral
st_dumppoints(pop.geom) p on true
) pops on pops.path = p
where
s < 360
group by
d.s
),
city_locs as (
select
d.s,
c.name,
st_setsrid(st_makeline(st_point(0,90,4326), p),4326) as geom,
deg,
case
/* label quadrant for qgis*/
when deg >= 337.5 and deg < 22.5 then 1
when deg >= 22.5 and deg < 67.5 then 2
when deg >= 67.5 and deg < 112.5 then 5
when deg >= 112.5 and deg < 157.5 then 8
when deg >= 157.5 and deg < 202.5 then 7
when deg >= 202.5 and deg < 247.5 then 6
when deg >= 247.5 and deg < 292.5 then 3
when deg >= 292.5 and deg < 337.5 then 0
end as lbl_box
from
cities c,
ellipsis d
join lateral
st_lineinterpolatepoint(
d.geom,
(round(c.lon,0)+180) / 360.0
) p on true
join lateral
degrees(
st_azimuth(
st_point(0,0,4326),
p
)
) deg on true
)
select
row_number() over()::int as oid, s, name,
('2022.11.20 '||(((24.0 * s) / 360.0||' hours')::interval)::time)::timestamp as ts,
geom, deg, lbl_box, cl
from (
select
s, geom, deg, lbl_box,
case
when lbl_box in (3,5,6,7,8) then name
else ''
end as name, 'cities' as cl
from
city_locs
union all
select
s, geom, null::numeric as deg, null::int as lbl_box,
null::varchar as name, 'pop_lines' as cl
from
pop_lines
) d
;