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

alphabetasoup
  • 8,718
  • 4
  • 38
  • 78
James
  • 305
  • 3
  • 7
  • For additional ideas see related threads at http://gis.stackexchange.com/questions/429 and http://gis.stackexchange.com/questions/18562. Conceivably the ideas presented in http://gis.stackexchange.com/questions/17788 could be helpful, too. I wonder, though, what you mean by "interior": these polygons are not well-defined, so at a minimum you need to supply information to indicate (a) which side of each polyline is considered "interior" and (b) how to cut them off near the poles. – whuber Apr 13 '16 at 16:30

3 Answers3

5

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.

AndreJ
  • 76,698
  • 5
  • 86
  • 162
3

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

MarMat
  • 181
  • 1
  • 2
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

Mr Purple
  • 1,451
  • 12
  • 25