Naive räumliche Interpolation in Postgres

import psycopg2
from multiprocessing import Pool

def route(modulo):
    connection = psycopg2.connect("service=osm")
    cursor = connection.cursor()
    query_string = "insert into routing.routes " \
        "with select_edge as ( " \
        "select routing_id, ST_X(ST_Transform(b.way, 4326)) || ',' || ST_Y(ST_Transform(b.way, 4326)) || ';' || " \
        "ST_X(ST_Transform(c.way, 4326)) || ',' || ST_Y(ST_Transform(c.way, 4326)) as route from " \
        "routing.edges as a, public.planet_osm_point as b, public.planet_osm_point as c " \
        "where b.osm_id=a.id1 and c.osm_id=a.id2 and routing_id % 10 = " + str(modulo) + " ) select routing_id, " \
        "(content::json->'routes'->0->>'distance')::numeric as distance, " \
        "(content::json->'routes'->0->>'duration')::numeric as duration from select_edge left join lateral " \
        "http_get('http://127.0.0.1:5000/route/v1/driving/' || route || '?steps=false&overview=simplified&annotations=false') as b " \
        "on true"

    cursor.execute(query_string,)
    connection.commit()
    cursor.close()
    connection.close()
    return True


if __name__ == '__main__':
    with Pool(10) as p:
        print(p.map(route, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]))

shp2pgsql -s 3035 STATISTIK_AUSTRIA_L001000_LAEA.shp statistik.regionalstatistische_rastereinheiten_001000_2016 | psql service=osm
shp2pgsql -s 3857 URAU_RG_100K_2020_3857_CITIES.shp eurostat.urban_audit_2020_cities | psql service=osm

with select_geom_ways as (
    select osm_id, way from eurostat.urban_audit_2020_cities as a, public.planet_osm_point as b 
    where cntr_code = 'AT' and place in ('town', 'city', 'village') and st_within(way, geom)
),

select_routes as (
    select id1, id2, duration  from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id1=c.osm_id
    union 
    select id2, id1, duration from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id2=c.osm_id
),

select_overall as (
    select id2 as osm_id, way, min(duration / 60) as duration from select_routes as a, public.planet_osm_point as b where a.id2=b.osm_id and id2 not in (
    select osm_id from select_geom_ways
)
group by 1,2
union
select osm_id, way, 0::numeric as duration from select_geom_ways

select_raster as (
    select id, st_transform(geom, 3857) as geom from statistik.regionalstatistische_rastereinheiten_001000_2016
),

select_osm_points as (
    select id, min(duration) as duration from select_raster as a,
    select_overall as b where ST_Within(b.way, a.geom)
    group by 1
)

get_neighbours as (
    select id, geom as geom, ST_LENGTH(ST_MakeLine(ST_Centroid(geom), way)) as distance, duration, osm_id from
    select_raster as a join select_overall as b on (ST_DWithin(a.geom, b.way, 15000))
    where id not in (select id from select_osm_points)
)

select id, avg(duration) as duration, percentile_cont(0.5) WITHIN group (order by duration) as duration_median 
from get_neighbours group by id
UNION
SELECT id, duration, duration FROM select_osm_points;

select a.id, geom, duration, duration_median from statistik.regionalstatistische_rastereinheiten_001000_2016 as a, 
spielwiese.osm_ort_duration as b where a.id=b.id

create table spielwiese.osm_ort_duration as
-- in OSM alle Punkte auswählen, die in die Urban Audit Klassifikation City fallen
with select_geom_ways as (
select osm_id, way from eurostat.urban_audit_2020_cities as a, public.planet_osm_point as b
where cntr_code = 'AT' and place in ('town', 'city', 'village') and st_within(way, geom)
),
-- alle routen von und zu diesen Punkten auswählen
select_routes as (
select id1, id2, duration  from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id1=c.osm_id
union
select id2, id1, duration from routing.edges as a, routing.routes as b, select_geom_ways as c where a.routing_id=b.routing_id and id2=c.osm_id
),
-- für jeden Punkt die minimale Fahrtdauer auswählen, die Dauer in Minuten umwandlen und die "City"-Punkte auf 0 setzen
select_overall as (
select id2 as osm_id, way, min(duration / 60) as duration from select_routes as a, public.planet_osm_point as b where a.id2=b.osm_id and id2 not in (
select osm_id from select_geom_ways
)
group by 1,2
union
select osm_id, way, 0::numeric as duration from select_geom_ways
),
-- den Raster de Statistik Austria auswählen und nach EPSG:3857 transformieren
select_raster as (
select id, st_transform(geom, 3857) as geom from statistik.regionalstatistische_rastereinheiten_001000_2016
),
-- die Raster auswählen, die schon einen Punkt beinhalten
select_osm_points as (
select id, min(duration) as duration from select_raster as a,
select_overall as b where ST_Within(b.way, a.geom)
group by 1
),
-- für alle Raster, die keine Punkte beinhalten, die "Nachbarn" berechnen
get_neighbours as (
select id, geom as geom, ST_LENGTH(ST_MakeLine(ST_Centroid(geom), way)) as distance, duration, osm_id from
select_raster as a join select_overall as b on (ST_DWithin(a.geom, b.way, 15000))
where id not in (select id from select_osm_points)
)
select id, avg(duration) as duration, percentile_cont(0.5) WITHIN group (order by duration) as duration_median
from get_neighbours group by id
UNION
SELECT id, duration, duration FROM select_osm_points;