-1

I have got a shapely polygon given in WGS84 coordinates formed like shown in the picture. As you can see, it intersects the date border.

enter image description here

When I get the polygon into a shapely polygon the coordinates are getting connected the wrong way and the polygon I get looks like this: enter image description here

How do I circumvent this issue?

Edit: This is how I load the polygon from a textfile, where it as given as POLYGON((-148.77203 44.959396,-147.72769 46.27871,- ...

   with open(footPrintFileName) as footprintFile:
        polygonString = footprintFile.read()
        filesFootPrintPolygon = shapely.wkt.loads(polygonString)

And this is how I the render it to a map

map = folium.Map(location=[51.077300, 10.205498], zoom_start=4)
folium.GeoJson(_geometry, style_function=lambda x: {'color': 'purple','weight': 1,'fillOpacity': 0.2},).add_to(map)
map.save(self.target_filepath)
webbrowser.open('file://' + self.target_filepath)

The rendering is not the problem tough. My target is to filter out some coordinates on the map based on that polygon, and if I do so only the polygons intersecting with the malicious polygon from the second image appear.

muilenheinz
  • 19
  • 3

1 Answers1

0

In the meantime, I found an completely different approach to my problem, but this is how far I came, in case anyone stumbles across this problem in the future: I just split the polygon into two polygons: Iterating over each pixel from the original polygon, at first all pixels are written into the first. As soon as the dateborder is crossed or one of the poles is "touched" I switch to the other and write the coordinates into there. When the location of the coordinates changes back, I change back to the first polygon. Now the coordinates at the edges of the map are missing and have to me inserted manually. This approach only works for a polygon with two subsections and has some other minor limitations, but it worked for my requirements.

This solution is not complete tough: For the approach I just explained, I need to know whether the original polygon "touches" the south or north pole. Currently I insert that information manually. If you have a solution to that, feel free to post it here.

# when the polygons get split into multiple parts exceeding the date border or top of the map, the coordinates
# representing the edges of the map are missing. This function determines which coordinate is missing at which place
# and fills it into the given polygon
def insertExtraPointsToPolygon(_polygon: Polygon):
    # check which position the new point has to be added at
    previousCoordinate = None
    maxDistance = 0
    storeCoordinates = None
    # append the first coordinate to the end, so the distance between those two will also be calculated
    # in the loop below
    polygonWithRepetition = _polygon
    polygonWithRepetition.append(_polygon[0])
    for coordinate in polygonWithRepetition:
        if previousCoordinate is not None:
            dist = geopy.distance.geodesic((coordinate[1], coordinate[0]), (previousCoordinate[1], previousCoordinate[0])).km
            if dist > maxDistance:
                maxDistance = dist
                storeCoordinates = previousCoordinate, coordinate
        previousCoordinate = coordinate
    index = _polygon.index(storeCoordinates[0]) + 1

    # check which point has to be added by checking which point has the lowest distance to the two "edge" coordinates
    # determined in the previous part of this function. This has to be done in coordinateDistance, since otherwise
    # the distance to e.g. [-180, 90] would be the same [180, 90], since those points are approximately on the same
    # spot on a globe ( but not on a map)
    min = float("inf")
    outputEdgePoint = None
    edgeCoordinates = [[-180, 90], [180, 90], [180, -90], [-180, -90]]
    for edgeCoordinate in edgeCoordinates:
        dist1 = abs(edgeCoordinate[1] - storeCoordinates[0][1]) + abs(edgeCoordinate[0] - storeCoordinates[0][0])
        dist2 = abs(edgeCoordinate[1] - storeCoordinates[1][1]) + abs(edgeCoordinate[0] - storeCoordinates[1][0])
        overallDistance = dist1 + dist2
        if overallDistance < min:
            min = overallDistance
            outputEdgePoint = edgeCoordinate

    _polygon.insert(index, outputEdgePoint)
    return _polygon


# most of the footprints overlap the date border. Shapely does not connect the coordinates correctly in this case, which
# leads to malicious footprint polygons. To circumvent this issue, spilt the polygon into two subsections
# @_polygon: Shapely polygon to be split
def correctPolygonShape(_polygon: Polygon):
    print("correct polygon shape...")
    coordinates = _polygon.exterior.coords
    previousCoord = None
    polygon1 = []
    polygon2 = []
    currentPolygon = polygon1
    max_xy = max(coordinates, key=lambda x: (x[1]))
    min_xy = min(coordinates, key=lambda x: (x[1]))

    # check whether the coordinates get close to the top or bottom border of the map
    # TODO: automated implementation of this step
    minOrMaxBorder = "min" # todo

    # loop all coordinates and assign them to their respective polygons
    for coordinate in coordinates:
        if previousCoord is not None:
            # get the intersection with the date border
            longitudeDistInDegrees = abs(coordinate[0] - previousCoord[0])
            if longitudeDistInDegrees > 300:
                currentPolygon = polygon1

            # get the "intersection" with the top of the map
            if (minOrMaxBorder == "max" and coordinate == max_xy) or (minOrMaxBorder == "min" and coordinate == min_xy):
                currentPolygon = polygon2

        currentPolygon.append(coordinate)
        previousCoord = coordinate

    polygon1 = insertExtraPointsToPolygon(polygon1)
    polygon2 = insertExtraPointsToPolygon(polygon2)

    poly1 = Polygon(polygon1)
    poly2 = Polygon(polygon2)
    footprint = MultiPolygon([poly1, poly2])
    # plot = Plot()
    # plot.plotGeometry(footprint)

    print("...done")
    return footprint
muilenheinz
  • 19
  • 3