10

I'm trying to generate polygons for satellite orbital swaths. So far I have a method to generate two lines which represent the edge of each swath in [lat,long]. Some of the swath's cross the international dateline and so wrap round:

swath wrap around

I was able to solve this with ogr2ogr -wrapdateline:

ogr2ogr -wrapdateline  -f "ESRI Shapefile" test.shp orbits.shp

Which does split the lines probably

I now want to be able to generate polygons on the interior of both lines. So for example in the case where one edge of the swath crosses the dateline a polygon fills in when it emerges on the other side, like:

fill

I need a method that is automated as I need to repeat the task a lot. Preferably in python as that's how i have generated the lines. Here are the two shapefiles containing the lines: wraparound ; datelinefixed

3

You can build a custom mercator projection centered approximately on the center of the swath. For example, use for swath 25:

+proj=merc +lon_0=-140 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

In this projection, the swath is not broken by the dateline. You can create the polygon from the line.

enter image description here

Then create a cut polygon between -179.95°E and 179.95°E in EPSG:4326:

Nr;WKT
1;POLYGON ((-179.95 89, 179.95 89, 179.95 -89, -179.95 -89, -179.95 89))

Reproject it to your custom CRS too, and subtract it from the swath polygon.

After reprojecting back to EPSG:4326, the swath is correctly divided by the dateline:

enter image description here

Continue with all swaths that cross the dateline.

| improve this answer | |
1

I would rewrite the swathe line generation process to start and finish in the same continuous longintudinal space. ie if a line started at 170° and finished at -170° I would rewrite the process to finish at 190° instead without wrapping at -180,180

Then you can make unbroken polygons between your lines.

Then use a clip process to split the polygons at the 180,-180 line and shift any parts which lie outside the -180,180 space by adding or subtracting 360° as appropriate.

Just get it all done before you save it with a particular projection/datum

| improve this answer | |
1

Thanks to @AndreJ for this idea, using Django GEOS API here is a simple solution that avoids needing to re-project anything:

1) Create a MultiPolygon that borders on the dateline:

from django.contrib.gis.geos.collections import MultiPolygon, LinearRing, Polygon
box1 = ((180.0, 89), (179.95, 89), (179.95, -89), (180.0, -89), (180.0, 89))
box2 = ((-180.0, 89), (-179.95, 89), (-179.95, -89), (-180.0, -89), (-180.0, 89))
poly1 = Polygon(LinearRing(box1))
poly2 = Polygon(LinearRing(box2))
poly = MultiPolygon(poly1, poly2)

2) If the offending geometry intersects, return the difference:

from django.contrib.gis.geos.geometry import GEOSGeometry
geometry = GEOSGeometry(WKT)  # WKT is your polygon in WKT string format
if geometry.intersects(poly):
    print("Geometry crosses dateline... splitting")
    geometry = geometry.difference(poly) # clip with dateline polygon

Result is shown as follows:

Before

After

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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