Seleccione el cuadro delimitador con postGIS?

36

Quiero crear una consulta para seleccionar todas las formas y sus nodos que existen dentro de un cuadro delimitador usando postGIS. El cuadro delimitador incluirá todos los detalles que el comando de ósmosis "--bounding-box" recuperará.

¿Hay alguna forma de hacer eso?

postgis openstreetmap vector osmosis extents uriel
fuente

Respuestas:

36

Para los documentos de osmosis, veo la opción de comando :

--bounding-box top=49.5138 left=10.9351 bottom=49.3866 right=11.201

para PostGIS puede usar ST_MakeEnvelope (izquierda, abajo, derecha, arriba, cuadrícula) para construir un cuadro delimitador, luego el &&operador del cuadro delimitador para encontrar dónde se intersecan los cuadros delimitadores:

SELECT *
FROM mytable
WHERE mytable.geom && ST_MakeEnvelope(10.9351, 49.3866, 11.201, 49.5138, 4326);

El SRID 4326 es para WGS84 Lat / Long, y solo se requiere para PostGIS 1.5; Se puede omitir para versiones posteriores.

Mike T
fuente
Gracias. la función ST_MakeEnvelope necesita un parámetro más, srid. No sé qué poner allí ... ¿alguna idea?
uriel
1
Parece que está utilizando PostGIS 1.5, que requiere ese parámetro. Creo que SRID se ignora, por lo que cualquier valor puede producir los mismos resultados. Si tiene datos de lat / long, generalmente use un SRID de 4326.
Mike T
1
La mayoría de las herramientas en estos días le permiten elegir el SRID para datos OSM cuando lo carga. El SRM de OSM predeterminado es 3857 (mercator esférico). El SRID predeterminado para la mayoría de los datos lat / lon es SRID 4326 (Lat / Lon AKA WGS84). Si carga los datos con SRID 3857, por ejemplo, tendrá que hacer una conversión de LAT / LON WGS84 a 3857: ST_Transform (ST_MakeEnvelope (LON1, LAT1, LON2, LAT2, 4326), 3857) Algunas herramientas (como imposm3) actualmente solo admite SRID 3857
Justin Swanhart
Tenga en cuenta que el operador && no transforma los SRID por usted. Asegúrese de que la envolvente que realice esté en el mismo SRID que la geometría de prueba, o bien transfórmela usted mismo. trac.osgeo.org/postgis/ticket/2320
Nelson
1
El operador && es computacionalmente más lento que ST_Intersects
caiohamamura
8

Creo que será algo como esto: el cuadro delimitador en PostGIS es creado por 

ST_GeomFromText('POLYGON((ulx uly, urx ury, llx llr, lrx lry, ulx uly))', <srid>)

La consulta utilizará ST_Intersection con una subconsulta.

SELECT bbox_nodes.id, bbox_nodes.tag, nodes_geom 
FROM (SELECT nodes.id, nodes.tag, 
   ST_Intersection(nodes.the_geom, 
      ST_GeomFromText('POLYGON((ulx uly, urx ury, llx llr, lrx lry, ulx uly))', <srid> )).geom AS nodes_geom
   FROM nodes 
   WHERE ST_Intersects(nodes.the_geom, 
      ST_GeomFromText('POLYGON((ulx uly, urx ury, llx llr, lrx lry, ulx uly))', <srid> )) AS bbox_nodes
WHERE ST_Dimension(bbox_nodes.nodes_geom)=0;

Más o menos tomé esto de las páginas de ayuda de PostGIS.
Una segunda consulta, en la tabla de formas, diseñada de manera similar a la anterior (pero con ST_Dimension () = 1) debería obtener las formas.

HTH, Micha

Micha
fuente
¡Hola Gracias! que diablos ¿Qué debo insertar en <srid>? y ".geom" (línea 4) parece no ser válido, ¿debería estar allí?
uriel
Lo siento, me perdí tu comentario de la semana pasada. La cuadrícula es el código del Sistema de referencia de coordenadas. es decir, 2039 para Israel. La adición .geom extrae la parte de Geometría de una "Colección de Geometría. Puede que tengas razón en que no se requiere aquí.
Micha
5

Hay un tema aquí que es similar a su pregunta aquí ...

ST_Intersection - (T) Devuelve una geometría que representa la porción compartida de geomA y geomB. La implementación de geografía hace una transformación a la geometría para hacer la intersección y luego transformarse nuevamente a WGS84.

1. También puede obtener información aquí sobre las funciones de construcción de geometría.

SELECT ST_AsText(ST_Intersection(
  ST_Buffer('POINT(0 0)', 2),
  ST_Buffer('POINT(3 0)', 2)
));

intersección

2.Otra información aquí sobre Intersección Intersección: PostGIS - ST_Intersects, ST_Intersection ...

SELECT b.the_geom As bgeom, p.the_geom As pgeom, 
        ST_Intersection(b.the_geom, p.the_geom) As intersect_bp
    FROM buildings b INNER JOIN parcels p ON ST_Intersection(b,p)
    WHERE ST_Overlaps(b.the_geom, p.the_geom)
    LIMIT 1;

intersección

Espero que te ayude...

Aragón
fuente
0

Este es un comentario sobre el código de @ Micha.

Los pares de coordenadas para el POLYGONdeben seguir un orden en el sentido de las agujas del reloj (o en sentido contrario): superior izquierda, superior derecha, inferior derecha, inferior izquierda, superior izquierda nuevamente.

Entonces, en sentido horario, la llamada a la función debería ser:

ST_GeomFromText('POLYGON((ulx uly, urx ury, lrx lry, llx llr, ulx uly))', <srid>)

O en sentido antihorario:

ST_GeomFromText('POLYGON((ulx uly, llx llr, lrx lry, urx ury, ulx uly))', <srid>)
Daishi
fuente