La unión espacial Postgres/PostGIS devuelve resultados nulos

Question

Tengo 1,3 millones de puntos y 15 mil polígonos.

Estoy tratando de hacer una unión espacial para saber cuántos puntos se cruzan con cada polígono. En última instancia, voy a sumar los valores de algunas columnas en los datos de puntos por cada polígono.

Así que los datos del punto son PCD4 y los del polígono son HCBG. He comprobado que ambas tablas tienen el mismo SRID. He probado esta consulta de unión espacial bastante sencilla, pero los resultados son un conjunto nulo que no tiene sentido.

select geoid, id
from public."PCD4"
join public."HCBG"
ON ST_Intersects(public."HCBG".geom, public."PCD4".geom)
group by id, public."PCD4".geom, geoid`

He probado a limitar la consulta a un solo polígono pero eso también se anula. He comprobado los datos en QGIS y definitivamente hay puntos en los polígonos y específicamente en el único polígono de prueba que probé.

Si intento hacer la unión espacial en QGIS, da error.

He comprobado la integridad de los polígonos. Son archivos de Census Tiger y han comprobado que están bien.

Así que estoy perdido. Estoy seguro de que el problema es algo que tiene que ver con la alineación de los datos, ya sea la proyección o algo así. Pero no estoy seguro de cómo solucionar el problema más allá de comprobar el SRID.

Answer 1

Cuanto más utilizo SQL, más me parece que es más sencillo dividir estas consultas en sus componentes, es decir, subconsultas, en lugar de intentar hacerlo todo de una vez.
Empecemos con una unión para etiquetar cada punto con el polígono en el que se encuentra. Estoy utilizando el método st_covers para esto, ya que tiene en cuenta los puntos que caen en el límite de un polígono. Asigno un alias a cada nombre de tabla para facilitar su referencia en la consulta. No estoy seguro de qué ID va con cada capa, así que me invento los nombres. Basta con sustituirlos por los correctos.

SELECT pcd4.fid,
       hcbg.gid,
       pcd4.geom
FROM public."PCD4" pcd4
JOIN public."HCBG" hcbg
  ON st_covers(hcbg.geom, pcd4.geom)

Esto selecciona el ID de cada punto, y asigna el ID del polígono que lo contiene. Esto se puede utilizar ahora como un WITH consulta para otro SELECT consulta.

WITH pcd4_poly AS (
  SELECT pcd4.fid,
         hcbg.gid
  FROM public."PCD4" pcd4
  JOIN public."HCBG" hcbg
    ON st_covers(hcbg.geom, pcd4.geom)
)
SELECT hcbg.gid,
       count(pcd4_poly.fid),
       hcbg.geom
FROM public."HCBG" hcbg
WHERE hcbg.gid = pcd4_poly.gid
GROUP BY hcbg.gid

Esta consulta selecciona el ID y la geometría de cada polígono. La dirección WITH ejecuta esa consulta select para seleccionar el ID de cada punto, etiquetado con el ID de cada polígono. Cuando se seleccionan estos juntos, se crea una selección de uno a muchos de los polígonos a los puntos. Agrupando el ID del polígono con un recuento del ID del punto se obtiene el resumen de registros de puntos por polígono que se está buscando.

Admito que esta puede no ser la forma más eficiente de hacerlo. Sin embargo, te da un proceso en el que es fácil probar cada componente para asegurarse de que funcionan por separado, antes de ponerlos juntos.

Para abordar lo que estaba en cuestión con su consulta, me doy cuenta de que incluye lo que creo que son los ID de ambas capas, además de la geometría de la capa de puntos en su GROUP BY cláusula. Dado que está haciendo una intersección, no hay necesidad de tener una cláusula de grupo ya que no está realizando ninguna operación de agregación en el SELECT consulta. La primera parte de su consulta sería suficiente. Es posible que haya tenido una NULL resultado simplemente porque intentaba mostrar 1,3 millones de puntos.

select geoid, id
from public."PCD4"
join public."HCBG"
ON ST_Intersects(public."HCBG".geom, public."PCD4".geom)