Bueno, ante todo quiero agradecer a Fernando Aguilar de Geoprop que me dejó compartir esto con ustedes. Lo que les voy a mostrar es como a partir de una tabla geográfica (postgis) de segmentos de calles vamos a armar una tabla con las intersecciones de esas calles. Seguramente ya se estarán haciendo la misma pregunta que yo le hice a Fernando ¿para qué?, o sea, porque insertar redundancia en el sistema, y la respuesta es muy simple, performance. Si bien no creo que nadie tenga este problema puntualmente es un buen ejemplo de como jugar con postgis.
Así que bueno, ahí vamos, lo primero que les muestro es las tablas que usaremos:
geoprop=# \d segmentos_calle;
Table "public.segmentos_calle"
Column | Type | Modifiers
------------------+----------------------+----------------------------------------------------------------------------
idsegmentoscalle | integer | not null default nextval('segmentos_calle_idsegmentoscalle_seq'::regclass)
idcalles | integer | not null
alturadesdeizq | smallint |
alturahastaizq | smallint |
alturadesdeder | smallint |
alturahastader | smallint |
anchoacera | real |
anchovereda | real |
cpa | character varying(7) |
the_geom | geometry |
Indexes:
"segmentos_calle_idsegmentocalle" PRIMARY KEY, btree (idsegmentoscalle)
"geomindex" gist (the_geom)
"segmentos_calle_callessegmento_calle" btree (idcalles)
"segmentos_calle_idcalles" btree (idcalles)
Check constraints:
"enforce_dims_the_geom" CHECK (ndims(the_geom) = 2)
"enforce_geotype_the_geom" CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL)
"enforce_srid_the_geom" CHECK (srid(the_geom) = 4326)
Foreign-key constraints:
"callessegmento_calle" FOREIGN KEY (idcalles) REFERENCES calles(idcalles) ON UPDATE CASCADE ON DELETE CASCADE
geoprop=# \d barrios;
Table "public.barrios"
Column | Type | Modifiers
-----------+-----------------------+-------------------------------------------------------
gid | integer | not null default nextval('barrios_gid_seq'::regclass)
id | bigint |
nombre | character varying(20) |
area | double precision |
perimetro | double precision |
the_geom | geometry |
Indexes:
"barrios_pkey" PRIMARY KEY, btree (gid)
"barrios_the_geom_gist" gist (the_geom)
Check constraints:
"enforce_dims_the_geom" CHECK (ndims(the_geom) = 2)
"enforce_geotype_the_geom" CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL)
"enforce_srid_the_geom" CHECK (srid(the_geom) = 4326)
Esto es, tenemos dos tablas, una con los segmentos de calles, y otra con los barrios. El objetivo es crear una tercer tabla con los ids de los segmentos que se interccionan en cada barrio, osea, idsegmento1, idsegment2 e idbarrio.
Antes que nada creamos la tabla:
CREATE TABLE intersecciones ( idsegmentoscalle1 integer NOT NULL, idsegmentoscalle2 integer NOT NULL, id_barrio bigint NOT NULL, CONSTRAINT "intersec_PK" PRIMARY KEY (idsegmentoscalle1, idsegmentoscalle2, id_barrio) ) WITH ( OIDS = TRUE ) ; ALTER TABLE intersecciones OWNER TO devgeoprop;
Bueno, ahora necesitamos pensar en lo que tenemos que conseguir.
1) Conseguir la lista de ides de segmentos que se interseccionan
2) conseguir el punto en donde lo hacen
3) conseguir el barrio en donde lo hacen
A por ello!
Vamos a crear dos funciones para evitar hacer selects anindados que al final complican la lectura:
1) getbarriofromgem(the_geom)
Nos devuelve el id de barrio de una determinada geometría
2) segmentosinter (idsegmento1, idsegmento2)
Nos devuelve la geometria de la intersección de dos calles.
CREATE OR REPLACE FUNCTION getbarriofromgem(geometry) returns bigint AS 'select id from barrios where ST_Contains(the_geom,$1);' LANGUAGE SQL IMMUTABLE RETURNS NULL ON NULL INPUT; CREATE or REPLACe FUNCTION segmentosinter(integer, integer) returns geometry AS 'select distinct Intersection(calle1.the_geom,calle2.the_geom) from segmentos_calle calle1, segmentos_calle calle2 where ST_Touches(calle1.the_geom, calle2.the_geom) and calle1.idsegmentoscalle = $2 and calle2.idsegmentoscalle = $1;' LANGUAGE SQL IMMUTABLE RETURNS NULL ON NULL INPUT;
Y ahora solo resta hacer nuestro simple select
insert into intersecciones (idsegmentoscalle1, idsegmentoscalle2, id_barrio) select distinct calle1.idsegmentoscalle, calle2.idsegmentoscalle, getbarriofromgem(segmentosinter(calle1.idsegmentoscalle, calle2.idsegmentoscalle)) from segmentos_calle calle1, segmentos_calle calle2 where calle1.the_geom && calle2.the_geom and calle1.idcalles <> calle2.idcalles;
Más simple imposible, esto se puede optimizar si se complica el select, metiendo las funciones dentro. Pero, es algo que se hace una sola vez para cargar la tabla, por lo tanto prefiero que se lea bien ;).
Ahora la tabla está cargadita :D, pero.. ¿qué pasas si se borra un segmento de calle?
Bueno, basta con añadir una restricción foraña:
ALTER TABLE intersecciones ADD CONSTRAINT "idsegmentos1_FK" FOREIGN KEY (idsegmentoscalle1) REFERENCES segmentos_calle (idsegmentoscalle) ON DELETE CASCADE; ALTER TABLE intersecciones ADD CONSTRAINT "idsegmentos2_FK" FOREIGN KEY (idsegmentoscalle2) REFERENCES segmentos_calle (idsegmentoscalle) ON DELETE CASCADE;
Ahora, ¿qué pasa si inserta un nuevo segmento?
Bueno, para eso vamos a crear un trigger:
create or replace function actualiza_segmentoscalle() returns trigger as $emp_stamp$ begin insert into intersecciones (idsegmentoscalle1, idsegmentoscalle2, id_barrio) select distinct calle1.idsegmentoscalle, calle2.idsegmentoscalle, getbarriofromgem(segmentosinter(calle1.idsegmentoscalle, calle2.idsegmentoscalle)) from segmentos_calle calle1, segmentos_calle calle2 where calle1.the_geom && calle2.the_geom and calle1.idcalles <> calle2.idcalles and calle1.idsegmentoscalle = NEW.idsegmentoscalle; return NEW; end; $emp_stamp$ LANGUAGE plpgsql; CREATE TRIGGER actualiza_segmentos_calle BEFORE INSERT ON segmentos_calle FOR EACH ROW EXECUTE PROCEDURE actualiza_segmentoscalle();
Con esto nos aseguramos que cuando se inserta un segmento nuevo se actualicen las intersecciones también.
Espero les sirva, peguenle una mirada al proyecto geoprop, quizás para cuando lean este artículo estará terminado.
Si necesitan un hosting de postgres con soporte de mapserver, geoserver y postgis les recomiendo Soluciones Root.
Saludos.
Tags: funciones postgis, hosting postgis, postgis, postgis ejemplo, postgis intersección, postgis sql, postgis triggers, postgres
Hola,
estoy intentando insertar una capa (fichero de extension .shp) con quantum gis a una base de dtaos postgres. Al insertar el fichero da el siguiente error:
Problemas al insertar objetos espaciales del archivo:
C:\Documents and Settings\user\Escritorio\pruebas\fichero_p.shp
La base de datos dio un error mientras ejecutaba esta SQL:
INSERT INTO “public”.”fichero_p” VALUES(0,’ 110000′, I’,’0′,’471.649′,NULL,NULL,NULL,’0′,… (cortado el resto de la SQL)
El error fue:
ERROR: new row for relation “fichero_p” violates check constraint “enforce_dims_the_geom”
Alguien me puede ayudar?.
Muchas gracias.
Laura, el problema es que está insertando una gemotría de 3 dimensiones en una tabla de dos dimensiones. Por eso viola el contraint, te recomiendo que uses el comando shp2psql. Te acabo de hacer un howto de como acerlo, puedes entrar en este post.
Espero te sirva.
Gracias por comentar.
Hola,
muchísimas gracias por contestar!!!.
Importe el fichero shp a la base de datos postgres como me indicaste pero el servidor de mapas geoserver no pilla la capa insertada.
Puede ser por la geometría de 3 dimensiones?.
El fichero shp que no puedo insertar con quantum gis porque tiene una geometría de 3 dimensiones me lo tiene que dar un chico, si le digo a este chico que me mande el fichero con una geometría de 2 dimensiones lo podre insertar con quantum gis?.
Muchas gracias.
Un saludo.
Laura.
Hola Laura, que el geoserver no pille la capa puede ser por muchas cosas, si te creo la tabla, y la miras y tiene la columna the_geom con algo dentro, algo debería dibujarte. Si no te dibuja nada se puede a que la proyección que sea diferente, entonces si te dibuja pero lo hace en cualquier lado y te da la impresión de que no funciona (creo que es el error más común). hay una tabla llamada spatail_ref_sys (creo) que tiene la lista de todos los SRID, si no puedes buscarlos por internet, la más común es la 4326 (WGS84), otra forma de ver si está todo bien es hacer un select y ver en que coordenadas te está poniendo el polígono, eso lo puedes hacer con la siguiente consulta:
Lo de quatum, no creo que te esté pasando el shape en tres dimensiones, lo que puede ser es que quatum esté teniendo un problema para exportar la capa a postgis, quizás un bug, podrías probar con otro programa por ejemplo uDig, aunque el shp2psql es lo mejor que hay ;).
Mi recomendacion es:
Averigua que proyección tiene tu shape, con quatum por ejemplo, y luego ejecuta el mismo comando pero agrega el parámetro -s
SELECT astext(the_geom) from mitabla;
Ahí te deberías dar cuenta si hay un problema de proyección como yo sospecho.
Gracias por tu comentario, una vez más espero ser de ayuda.
MN
Hola,
muchas gracias por reponderme.
Sigo con problemas con los ficheros shp.
Pase el fichero shp a un fichero sql asi: Shp2pgsql -s 23030 C:\prueba\carpeta\fichero.shp public.prueba > C:\prueba\carpeta\ficheroSQL.sql pero el goserver no lo pilla.
Despues hice lo q me dijiste de SELECT astext(the_geom) from mitabla pero me sale un churro de cosas como MULTILINESTRING(………..) y no se las coordenadas que esta pintando.
Como puedo saber la proyección que tiene el fichero?.La proyección que estoy utilizando es 23030 que es la de España.
Muchas gracias.
Laura.
Hola Laura, prueba usar ogr. Con el siguiente comando:
ogr2ogr -overwrite -nlt NONE -lco PRECISION=NO -f PostgreSQL -t_srs “EPSG:4326” -s_srs “EPSG:22194” “PG:dbname=granchaco user=pgsql” TUARCHIVO.shp
-t_srs es la proyección destino, (siempre es bueno tener todo en wgs84).
-s_srs es la proyección de shape, que podrías sacarla con un ogrinfo, o bien obviar ese parámetro.
Otra cosa, te mando un correo con mi correo, pásame por ahí la salida del select que haces, a ver si a simple vista veo la proyección.
Avísame si tienes éxito.
MN
Hola,
muchísinas gracias!!!, ya me funciona, era un tema de proyecciones.
La verdad es que estoy muy pez en temas de cartografía. Te voy hacer
otra pregunta haber si me puedes ayudar, l. Tengo en un mapa dos
capas de lineas, tengo que diferencias las capas, una ponerla en color
verde y otra en color rojo, sabes como puedo cambiar el estilo de la
linea de cada capa con el geoserver?.
Muchas gracias. De todos los comentarios que puse eres la única persona que me contesto.
Un saludo.
Laura
Muchas gracias!!!, ya conseguir cambiar el color.
Repito muchas gracias por todo.
Un saludo.
Laura.
Me alegro mucho Laura, devuelvele el favor a otro que esté en apuros 😉
Un beso.
MN
Hola,
te voy hacer otra pregunto, el geoserver al interpretar los ficheros los interpreta mal, en lugar de mostrar lineas muestas manchas en el mapa, nada que se parezca a lineas. Sabes que puede ser?.
Lo más raro es que en otro ordenador me funciona correctamente. Comprobe los datos con la sentencia
SELECT astext(the_geom) from mitabla; y son iguales en los dos ordenadores.
Muchas gracias.
Laura.
Hola Laura, es posible que esto tenga que ver con el estilo del geoserver, quizás tiene un estilo de punto o de polígono en lugar de tener un estilo de línea.
Espero te sirva.
MN
Hola,
Antes de nada muchas gracias.
El estilo que tengo en el geoserver es line, por ejemplo si pongo point no funciona. El geoserver lo tengo en un tomcat y cuando quiero que muestre la capa en la consola del tomcat sale:
Posible uso de la proyeccion “Tranverse_Mercator ” fuera de sua area de validez.
La latitud esta fuera de los limites permitidos.
Tienes idea que puede estar pasando?.
Muchisimas gracias.
Laura.
Hola Laura,
mmmm puede ser un error al importar los shapes con otra proyección. Si podes mandame el shape al mi correo y le doy una mirada con el comando posta para que salga. Por los datos que me das acá no se me ocurre que puede ser. Adjúntame también una captura de tu imagen de geoserver.
Lamento no poder ser de más ayuda por acá :(.
MN
Hola, estoy intentando cargar capas postgis en Geoserver, no he tenido problemas con capas shp pero creo un almacen postgis, publico las capas y les asigno un estilo y cuando intento cargar el servicio WMS no me las carga. Como se pueden cargar?. Al crear el almacen le doy todos los datos: el servidor, la base de datos, el puerto, el usuario y la clave y.. nada que no funciona.
Hola Juan, quizás lo que te está faltando es reiniciar el geoserver. Si te muestra el servicio listado lo más común es que sea un problema de proyecciones, que está mostrando realmente pero en otro lugar de la tierra, puedes probar especificando la proyección correcta en el layer.
Saludos, gracias por el comentario.
MN
Hola Matías, tengo el mismo problema que Laura, importé desde Qgis una capa shp a una BD postgis , luego le agregué 2 columnas ; ahora la quiero ver desde mapserver , y el resultado es que se queda cargando eternamente.
Las proyecciones son: (esto lo consulto desde qgis)
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
y en mapserver tengo configurado:
“proj=longlat”
“ellps=WGS84″
“datum=WGS84″
“no_defs”
Lo cierto es que no me da error y queda cargando.
Muchas Gracias , muy buenas explicaciones.
Hola Jorge, prueba con ogr2ogr debería funcionarte igual que le funcionó a Laura.
Cualquier cosa me avisas.
MN
Hola Matias, en caso de tener cargadas las calles en una tabla “Calles” y no querer hacer la tabla relacional de intersecciones, como seria el sql para calcular entre que calles se encuenta la ubicacion?
Los campos en la tabla son label, type y geom.
Desde ya mil gracias!
Pablo
Hola Pablo, ¿cual sería el where? o sea, ¿cual es el dato que vos tenés y en base a eso qué datos estás queriendo obtener? Decime eso y te armo un post 😉
Saludos y gracias por el comentario
Gracias por el post nuevo. El dato que yo tengo es un punto que deberia estar posicionado sobre una calle (o sea un punto sobre la linea) y quiero poder calcular en base a ese dato, entre que calles se encuentra la calle que contiene al punto original.
Saludos y gracias nuevamente.
Pablo
Hola Pablo, podes usar la función touch si es exacta, o bien la función touch combinada con la función buffer si es inexacta, por ejemplo:
SELECT id_calles from calles where calles.the_geom && astext(‘POINT (XXXX,XXXX)’); o si fuere inexacta:
SELECT id_calles from calles where calles.the_geom && buffer(astext(‘POINT (XXXX,XXXX)’),MARGEN);
La unidad de margen está en la misma unidad que tu sistema de proyecciones, usas WGS84 (SRID 4326) serán degrees, o si usas una proyección métrica serán metros.
Espero que te sirva, si tienés más dudas mándame un correo con tu schema y te hago un post.
Saludos.
Y en el caso en que hay calles que se llaman de un lado de la interseccion de una manera y del otro lado de la interseccion de otra, como calcular con sql que la continuacion del otro lado de la calle por mas que cambia el nombre, no es interseccion sino una continuacion con diferente nombre?
La verdad el blog me ha salvado mas de una vez.
Saludos y gracias de ante mando
Creo que la consulta que te dí no le importa el nombre de las calles, siempre y cuando estén en ids de calles distintas, es decir, si la calle A y B (B es continuación de A) tienen distintos IDs de calle, si no lo que tendrías que hacer es seleccionar los nombres de las calles en lugar de los IDs y luego joinear para buscar los IDs si es que los necesitas. ¿Me explico?
MN
Pablo lee este post, seguro que es lo que estás buscando.
Saludos y gracias de nuevo.
MN
hola admin, te felicito por las respuestas che, sos muy groso y buena persona respondiendo a todo. un saludo desde bsas, abrazo!
yo estoy lidiando con mapserver-openlayers-mapfish-ext-geoext, y realmente es medio un dolor de huevos todo esto! pero va saliendo.
estoy quedado en un buscador que me busque x direccion o unidad minima (parcela) y que me lleve al lat long de la parcela me la centre y me lleva a un nivel de zoom cercano con la parcela seleccionada.
puff
busco y busco y no encuentro che! si tenes alguna punta, muy agradecido!
estamos podiendo obtener un bbox de la unidad minima usando esto:
select smp, ST_box2d(the_geom) from parcelas where smp=’20-028-012′ ;
(nuestra unidad minima es seccion-manzana-parcela o sea s-m-p)
esto nos da el bbox (n,n,n,n) de la geometria, pero nos cuesta un poco estooooo abrazoooooooo 🙂
como catso obtenemos el lat long del centro de ese bbox? es necesario o nos estamos volviendo locos al dopes?
ja buenoooooooo amigazo,
un gran saludo, muy muy bueno tu blog!!!!!!
vladi
Hola Vladimir, mira hay una función que se llama centroid que te podría servir para sacar el punto del medio del bbox que seleccionaste, espero que te sirva.
Un saludo desde Corrientes.
MN
[…] a tunear un poco las consultas, sobre todo porque si no lo hacen las consultas que se usan en este post, serán muy lentas. Primero, creamos un índice en la […]
[…] estoy medio apurado así que voy al grano. Imaginen por ejemplo la tabla intersecciones de este post, pues anda todo muy bonito, pero tenemos el problema que se duplican las tuplas, por ejemplo: […]
hola buenos como estan soy novato en el tema y me gustaria saber si alguien sabe como sacar el punto central o centro de una geometry en postgres Gracias