PostGIS: interpretación incorrecta de un polígono que corta el meridiano 180

Question

Intento averiguar todos los geopuntos que intersecan el polígono establecido como parámetro.

El problema es cuando paso polígono que abarca aproximadamente una zona del estrecho de Bering (cerca de 180 longitud):

Así que uso la consulta:

SELECT ST_AsText(l.geo_point)
FROM "lightnings" "l"
WHERE (ST_Intersects(ST_GeomFromText('Polygon((132.0 40.0, -148.0 40.0, -148.0 -8.0, 132.0 -8.0, 132.0 40.0) )', 4326), geo_point));

Como puede ver, los vértices están colocados en el orden correcto, en el sentido de las agujas del reloj, desde el Noroeste. Pero el resultado cubre un área exterior e incluye todo otro mundo.

Por ejemplo, en resultado:

POINT(75.5637 40.0434)

El problema no toca el meridiano 0.

Prueba elemental:

SELECT ST_Area(ST_GeomFromText('Polygon((0.0 60.0, 10.0 60.0, 10.0 40.0, 0.0 40.0, 0.0 60.0) )', 4326))
UNION ALL
SELECT ST_Area(ST_GeomFromText('Polygon((-5.0 60.0, 5.0 60.0, 5.0 40.0, -5.0 40.0, -5.0 60.0) )', 4326))
UNION ALL
SELECT ST_Area(ST_GeomFromText('Polygon((175.0 60.0, -175.0 60.0, -175.0 40.0, 175.0 40.0, 175.0 60.0) )', 4326))

da el resultado:

200
200
7000

¿Hay algún truco sencillo para forzar que PostGIS me entienda? No me gusta la idea de dividir el polígono...

Answer 1

PostGIS tiene una función ST_Shift_Longitude(geom) para devolver geometrías Lon/Lat como 0-360 grados en lugar de -180 - 180.

Envuelve tus polígonos en esto, y los cálculos de área resultantes deberían ser correctos.

Descripción ici

BTW, supongo que te das cuenta de que tu cálculo de ST_Area está dando grados cuadrados. No es realmente útil... Si utiliza el nuevo tipo de datos de Geografía, (ST_GeographyFrom Text(...) ) entonces el cálculo del Área será en metros.

Adición:

Esto es lo que obtengo

geodata=# SELECT ST_Area(
   ST_GeomFromEWKT('SRID=4326;POLYGON((
   -179.5 0, 179.5 0, 179.5 1, -179.5 1, -179.5 0))')
   );

 st_area 
---------
     359
(1 row)

geodata=# SELECT ST_Area(
   ST_SHIFT_Longitude(
   ST_GeomFromEWKT('SRID=4326;POLYGON((
   -179.5 0, 179.5 0, 179.5 1, -179.5 1, -179.5 0))')));

 st_area 
---------
       1
(1 row)

Answer 2

Hice y respondí a una pregunta similar aquí: Círculo geométrico mínimo que cruza el meridiano 180º.

Lo que puedes hacer (aunque esto puede ser aún más feo que dividir el polígono) es utilizar valores superiores a 180/-180. Por ejemplo, en lugar de tener polygon((179, 0), (-179, 0)), puede utilizar polygon((179, 0), (181, 0)). La pregunta a la que he enlazado tiene un script en python que desplaza polígonos enteros de esta manera. Podrías adaptarlo para desplazar vértices dentro de tu polígono.