Hace unos años bastante extenso sobre intersecciones que se llamaba un poco de postigs, en este caso voy a intentar hacer lo mismo pero usando OSM, espero lo disfruten, a por ello!!!.
Primero creamos la tabla de intersecciones:
CREATE TABLE interseccion (
ID SERIAL PRIMARY KEY,
calle1 TEXT NOT NULL,
calle2 TEXT NOT NULL
);
SELECT AddGeometryColumn ('public','interseccion','geom',4326,'POINT',2);
Luego toca llenar la tabla…
insert into interseccion (geom,calle1,calle2) select distinct st_centroid(ST_Intersection(calle1.the_geom,calle2.the_geom)), calle1.name, calle2.name from ways calle1, ways calle2 where st_intersects(calle1.the_geom,calle2.the_geom) and calle1.name <> calle2.name and calle1.name <> '' and not calle1.name is null and calle2.name <>'' and not calle2.name is null;
Ahora borramos los duplicados ya que calle1 se intercepta con calle2, pero también calle2 con calle1, entonces borramos los que tienen el mismo punto:
DELETE FROM interseccion WHERE id IN (SELECT id FROM (SELECT id, ROW_NUMBER() OVER (partition BY geom ORDER BY id) AS rnum FROM interseccion) t WHERE t.rnum > 1);
Ahora creamos el índice y las columnas que faltan:
alter table interseccion add COLUMN ciudad text; alter table interseccion add COLUMN provincia text; alter table interseccion add COLUMN pais text; create index geom on interseccion using gist (geom);
y rellenamos con los datos correspondientes:
update interseccion set ciudad = (select name from planet_osm_polygon where way && st_transform(interseccion.geom, 900913) and boundary = 'administrative' and admin_level = '6'); update interseccion set provincia = (select name from planet_osm_polygon where way && st_transform(interseccion.geom, 900913) and boundary = 'administrative' and admin_level = '4'); update interseccion set pais = (select name from planet_osm_polygon where way && st_transform(interseccion.geom, 900913) and boundary = 'administrative' and admin_level = ‘2’);
para buscar usando el índice la consulta tendría que ser esta:
select *, st_astext(geom) from interseccion where ST_distance(geom, st_geomfromtext('point(-73.6388433333333 4.14556166666667)',4326))<0.0002 order by ST_distance(geom, st_geomfromtext('point(-73.6388433333333 4.14556166666667)',4326)) limit 1;
Y listo con eso lo tenemos, espero les sirva.
Saludos.
Tags: calles, intersecciones, nominatim, open street map, osm, postgis, reverse geocoding