Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I try to find out all geopoints which intersect the polygon set as a parameter.

The problem is when I pass polygon that roughly covers an area of Bering Strait (nearby 180 longitude): enter image description here

So I use the query:

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));

As you can see, vertexes are set in a correct order, clockwise, from the North-West. But the result covers outside area and including whole other world.

For example, in result:

POINT(75.5637 40.0434)

The problem doesn't touch the 0th meridian.

Elementary test:

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))

gives the result:

200
200
7000

Is there any simple trick to force PostGIS understanding me? I don't like an idea to divide the polygon...

share|improve this question

2 Answers 2

PostGIS has a built in function ST_Shift_Longitude(geom) to return Lon/Lat geometries as 0-360 degrees instead of -180 - 180.

Wrap your polygons in this, and the resulting area calculations should be OK.

Description here

BTW, I guess you realize that your ST_Area calculation is giving square degrees. Not really useful... If you use the new Geography data type, (ST_GeographyFrom Text(...) ) then the Area calculation will be in meters.

Addition:

Here's what I get

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)
share|improve this answer
    
Thanks! Of course I see what ST_Area means in this projection :) That's just for obviousness of PostGIS logic. –  Alexander Palamarchuk Jul 20 '12 at 19:25
    
Good function. But the result of test is almost the same, but wrong for the other meridian: 200, 700, 200 –  Alexander Palamarchuk Jul 20 '12 at 19:34
    
I've added a sample to the Answer above. Using ST_Shift_Longitude() returns the correct area for a polygon crossing the -180 deg boundary. –  Micha Jul 21 '12 at 15:39
    
Yes, you're right, but this situation is human-controlled. I gave to work with polygons generated by Google Maps - on the whole world. So ST_SHIFT_Longitude helps with the 180th meridian but plunge into the same troubles with the 0th meridian. Just try SELECT ST_Area(ST_Shift_Longitude(ST_GeomFromText('Polygon((-5.0 60.0, 5.0 60.0, 5.0 40.0, -5.0 40.0, -5.0 60.0))', 4326))) –  Alexander Palamarchuk Jul 21 '12 at 17:21
    
The only way I can think to work this out is make a line feature at the '180' meridian. Now create a function to test for each polygon if it intersects that line. If yes, then use the ST_Shift_longitude, otherwise the regular ST_Area. HTH –  Micha Jul 21 '12 at 19:16

I asked and answered a similar question here: The minimum bounding circle of geometry that crosses the 180th meridian

What you can do (though this may be even more ugly than dividing the polygon) is use values over 180/-180. E.g. rather than having polygon((179, 0), (-179, 0)), you can use polygon((179, 0), (181, 0)). The question I linked to has a python script that shifts whole polygons like this. You could adapt it to shift vertexes within your polygon.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.