The STs* of PostGIS
* saints?* spatial types?
preso/
balticgit-2022/
where's that bus at?
(supposedly)(data by the Estonian Road Administration)
length: ~ 3.5 km
Based on current time calculate fraction of time between trip start and end and then
st_lineinterpolatepoint
the same amount
length: ~ 3.5 km, total time: 2 minutes
(but with intermediate stops this bus will never be on time)
A to B: 1.8 km, B to A: 1.7 km
A to B: 30 secs, B to A: 90 secs
(add impedance to traveltime at stops)
Adding acceleration/deceleration to transit location calculation.
(see here for details)
(see here for details)
accelerate / decelarate: 8 secs, stop 3 + 3 secs
Red: interpolate whole trip. Blue: interpolate using stops. Green: interpolate using stops and impeded time
example dashboard and
API
the projection tango
madagascatorI think I've earned my geonerd points today.
— ⚛⚛⚛Iván Sánchez⚛⚛⚛ (@RealIvanSanchez) April 29, 2022
(re: https://t.co/84PGLpHLDZ ) pic.twitter.com/MZGwI76zSG
select
*
from (
select
gid, (st_dumppoints(geom)).*
from
ne_10m_admin_0_countries ) f
;
st_dumppoints for all vertices
select
st_project(
st_setsrid(st_point(-30,0), 4326)::geography,
st_distance(st_setsrid(st_point(0,90),4326)::geography, geom::geography),
st_azimuth(st_setsrid(st_point(0,90),4326)::geography, geom::geography) - radians(125)
)::geometry(point, 4326) as geom,
row_number() over()::int as oid, gid, path
from (
select
gid, (st_dumppoints(geom)).*
from
ne_10m_admin_0_countries ) f
where
st_y(geom) >-88
;
st_project
a new point (
geography
) with the old distance and azimuth (minus some whatever constant)
Based on NaturalEarth
country boundaries, and proccessed in PostGIS using
this SQL
select
*
from (
select gid, (st_dumppoints(geom)).*
from ne_10m_admin_0_countries
) f
;
st_dumppoints for all vertices
select
st_buffer(
st_setsrid(st_point(0,0),4326)::geography,
st_distance(st_setsrid(st_point(0,90),4326)::geography, geom::geography) / pi()
)::geometry(polygon, 4326)
from (
select gid, (st_dumppoints(geom)).*
from ne_10m_admin_0_countries
) f
;
st_project
will not be of use here, let's pretend we're on a flat surface so
st_buffer instead
select
st_rotate(
st_exteriorring(
st_buffer(
st_setsrid(st_point(0,0),4326)::geography,
st_distance(st_setsrid(st_point(0,90),4326)::geography, geom::geography) / pi()
)::geometry(polygon, 4326)
),
radians(90)
)
from (
select gid, (st_dumppoints(geom)).*
from ne_10m_admin_0_countries
) f
;
Take the resulting st_exteriorring
and st_rotate it 90 degrees
counter-clockwise
select
st_scale(
st_rotate(
st_exteriorring(
st_buffer(
st_setsrid(st_point(0,0),4326)::geography,
st_distance(st_setsrid(st_point(0,90),4326)::geography, geom::geography) / pi()
)::geometry(polygon, 4326)
),
radians(90)
), 1.70, 1.70
)
from (
select gid, (st_dumppoints(geom)).*
from ne_10m_admin_0_countries
) f
where st_y(geom) >-85 and st_y(geom) < 85
;
select
st_lineinterpolatepoint(
st_scale(
st_rotate(
st_exteriorring(
st_buffer(
st_setsrid(st_point(0,0),4326)::geography,
st_distance(st_setsrid(st_point(0,90),4326)::geography, geom::geography) / pi()
)::geometry(polygon, 4326)
),
radians(90)
), 1.70, 1.70
),
st_azimuth(st_setsrid(st_point(0, 90),4326)::geography, geom::geography) / (2*pi())
) as geom, gid, path
from (
select gid, (st_dumppoints(geom)).*
from ne_10m_admin_0_countries
) f
where st_y(geom) >-85 and st_y(geom) < 85
;
And finally st_lineinterpolatepoint
from the start of the exterior line to the fraction defined by
st_azimuth between the
North Pole and the original vertice.
Based on NaturalEarth
country boundaries, and proccessed in PostGIS using
this SQL
#30DayMapChallenge
(on GitHub and Twitter)(data courtesy of Estonian Environment Agency)
Unknown Pleasures by Joy Division
(and air temperature forecast)
image: beyondvinyl.co.uk
select
st_force3d(
st_point(
x,
y + temperature * 1000
)
)
from (
select temperature, valid_from, valid_to, x, y from phenomen
) f
;
create a point (st_point)
shifting its y coordinate to the north by coefficent of temperature, and st_force3d.
select
st_setsrid(
st_translate(
st_force3d(
st_point(
x,
y + temperature * 1000
)
), 0, 0, temperature
), 3301
)
from (
select temperature, valid_from, valid_to, x, y from phenomen
) f
;
select
st_makeline(
array_agg(
st_setsrid(
st_translate(
st_force3d(
st_point(
x,
y + temperature * 1000
)
), 0, 0, temperature
), 3301
) order by x
)
)
from (
select temperature, valid_from, valid_to, x, y from phenomen
) f
group by y, valid_from, valid_to;
st_makeline of all points on the same y
select
st_chaikinsmoothing(
st_makeline(
array_agg(
st_setsrid(
st_translate(
st_force3d(
st_point(
x,
y + temperature * 1000
)
), 0, 0, temperature
), 3301
) order by x
)
), 5
) as geom, y, valid_from, valid_to
from (
select temperature, valid_from, valid_to, x, y from phenomen
) f
group by y, valid_from, valid_to;
apply st_chaikinsmoothing on the lines
(sure, why not but for that we'll need to ...)
with segments as (
select
oid, st_dumppoints(geom) as pt
from
chaikinlines_temperature
)
select
*
from
segments
;
st_dumppoints
for all points of st_chaikinsmoothing-ed
lines.
with segments as (
select
st_makeline(
lag((pt).geom, 1, null) over(
partition by y, valid_from order by valid_from, y, (pt).path
),
(pt).geom
) as geom, valid_from, valid_to
from (
select
oid, st_dumppoints(geom) as pt
from
chaikinlines_temperature
) as dumps
)
select
*
from
segments
where
geom is not null;
st_makeline
of two consequtive points on the same y
with segments as (
select
st_makeline(
lag((pt).geom, 1, null) over(
partition by y, valid_from order by valid_from, y, (pt).path
),
(pt).geom
) as geom, valid_from, valid_to
from (
select
oid, st_dumppoints(geom) as pt
from
chaikinlines_temperature
) as dumps
)
select
geom, st_z(st_lineinterpolatepoint(geom, 0.5)) as temperature, valid_from, valid_to
from
segments
where
geom is not null;
take the midpoint of the segment with
st_lineinterpolatepoint
and its (interpolated)
st_z as temperature value
Polyps
(and wind speed + direction forecast)
- eccentricity for wind speed
- turn angle for wind direction
select
*
from (
select
st_transform(
st_project(
grid::geography, coalesce(wind_speed_mps,0)*500, radians(coalesce(wind_direction_deg, 0))-pi()
)::geometry(point, 4326),
3301) as geom, wind_speed_mps, wind_direction_deg, valid_from, valid_to
from
phenomen
) f;
shift grid point off its location by a factor of speed in the wind direction using
st_project and then
st_transform to
epsg:3301
select
st_ellipse(
st_x(geom), st_y(geom),
2500, 2500 + (coalesce(wind_speed_mps,0)*1000), 0, 16
)
from (
select
st_transform(
st_project(
grid::geography, coalesce(wind_speed_mps,0)*500, radians(coalesce(wind_direction_deg, 0))-pi()
)::geometry(point, 4326),
3301) as geom, wind_speed_mps, wind_direction_deg, valid_from, valid_to
from
phenomen
) f;
create a st_ellipse
in the previously shifted position with the longer axis as a factor of wind speed
select
st_rotate(
st_setsrid(
st_ellipse(
st_x(geom), st_y(geom),
2500, 2500 + (coalesce(wind_speed_mps,0)*1000), 0, 16
),
3301
),
-1*(radians(wind_direction_deg) + pi()),
geom
) as geom, wind_speed_mps, wind_direction_deg, valid_from, valid_to
from (
select
st_transform(
st_project(
grid::geography, coalesce(wind_speed_mps,0)*500, radians(coalesce(wind_direction_deg, 0))-pi()
)::geometry(point, 4326),
3301) as geom, wind_speed_mps, wind_direction_deg, valid_from, valid_to
from
phenomen
) f;
st_rotate
the ellipse in the desired direction anchoring on the shifted position.
(me going full-on estonian)
Polygonal expansion
to fill empty spaceReduce the number of vertices for coastline administrative units Zipcode areas (follow the next section on the right)
st_dumprings and
st_exteriorring on polygons to
get boundaries as linestrings
or ...
left is missing or right is missing
to make them a wee shorter at both ends
to extract the vertices ..
st_voronoipolygons
st_intersection with an
area-of-interest will give only those parts of Voronoi cells we're interested in
A closeup of the Voronoi cells colored by the polygon they belong to
st_union with the original self, and
st_difference with all other originals
A closeup of the final expansion product
from points to areas
(my little adventure with creating zip code ares)
Tartu: address point voronois colored by zipcode value
Tartu: roads, rivers, railways etc. lines
Tartu: road, river, railway etc. lines and overlaid addresspoints colored by zip
Tartu: road, river, railway etc. lines smashed into polygons (quartiers) and overlaid addresspoints colored by zip
All quartiers with only one zip
All cadastral parcels with only one zip
.. and merging these two together gives us already a pretty good coverage
Divide the empty space within quartiers that have zips assigned
Deal with address-pointless quartiers: longest bordering neighbor wins!
Zipcode areas (colored by zipcode value) with address point voronoi boundaries overlayed (red lines) for city of Tartu
Zipcode areas for the whole country
https://tkardi.ee/writeup/post/2018/07/22/from-points-to-polygons/
a shoutout to @tjukanov
https://media.ccc.de/v/bucharest-189-3-6-million-points-to-polygons-lessons-learned-while-generating-voting-districts-with-qgis-postgis-and-openjump-
as well useful in real life
tonis.kardi@gmail.com
https://tkardi.ee