-
-
Save PawaritL/ec7136c0b718ca65db6df1c33fd1bb11 to your computer and use it in GitHub Desktop.
| geojson = { | |
| "type": "Polygon", | |
| "coordinates": [ | |
| [ | |
| [179.0, 0.0], | |
| [-179.0, 0.0], | |
| [-179.0, 1.0], | |
| [179.0, 1.0], | |
| [179.0, 0.0] | |
| ] | |
| ] | |
| } |
| from typing import Union, List, Generator | |
| from shapely.geometry import mapping, Polygon, GeometryCollection | |
| from shapely import affinity | |
| def check_crossing(lon1: float, lon2: float, validate: bool = False, dlon_threshold: float = 180.0): | |
| """ | |
| Assuming a minimum travel distance between two provided longitude coordinates, | |
| checks if the 180th meridian (antimeridian) is crossed. | |
| """ | |
| if validate and any([abs(x) > 180.0 for x in [lon1, lon2]]): | |
| raise ValueError("longitudes must be in degrees [-180.0, 180.0]") | |
| return abs(lon2 - lon1) > dlon_threshold | |
| def translate_polygons(geometry_collection: GeometryCollection, | |
| output_format: str = "geojson") -> Generator[ | |
| Union[List[dict], List[Polygon]], None, None | |
| ]: | |
| for polygon in geometry_collection: | |
| (minx, _, maxx, _) = polygon.bounds | |
| if minx < -180: geo_polygon = affinity.translate(polygon, xoff = 360) | |
| elif maxx > 180: geo_polygon = affinity.translate(polygon, xoff = -360) | |
| else: geo_polygon = polygon | |
| yield_geojson = output_format == "geojson" | |
| yield json.dumps(mapping(geo_polygon)) if (yield_geojson) else geo_polygon |
| import math | |
| import copy | |
| import json | |
| from typing import Union, List | |
| from shapely.geometry import Polygon, LineString, GeometryCollection | |
| from shapely.ops import split | |
| # https://gist.github.com/PawaritL/ec7136c0b718ca65db6df1c33fd1bb11 | |
| from geopolygon_utils import check_crossing, translate_polygons | |
| def split_polygon(geojson: dict, output_format: str = "geojson", validate: bool = False) -> Union[ | |
| List[dict], List[Polygon], GeometryCollection | |
| ]: | |
| """ | |
| Given a GeoJSON representation of a Polygon, returns a collection of | |
| 'antimeridian-safe' constituent polygons split at the 180th meridian, | |
| ensuring compliance with GeoJSON standards (https://tools.ietf.org/html/rfc7946#section-3.1.9) | |
| Assumptions: | |
| - Any two consecutive points with over 180 degrees difference in | |
| longitude are assumed to cross the antimeridian | |
| - The polygon spans less than 360 degrees in longitude (i.e. does not wrap around the globe) | |
| - However, the polygon may cross the antimeridian on multiple occasions | |
| Parameters: | |
| geojson (dict): GeoJSON of input polygon to be split. For example: | |
| { | |
| "type": "Polygon", | |
| "coordinates": [ | |
| [ | |
| [179.0, 0.0], [-179.0, 0.0], [-179.0, 1.0], | |
| [179.0, 1.0], [179.0, 0.0] | |
| ] | |
| ] | |
| } | |
| output_format (str): Available options: "geojson", "polygons", "geometrycollection" | |
| If "geometrycollection" returns a Shapely GeometryCollection. | |
| Otherwise, returns a list of either GeoJSONs or Shapely Polygons | |
| validate (bool): Checks if all longitudes are within [-180.0, 180.0] | |
| Returns: | |
| List[dict]/List[Polygon]/GeometryCollection: antimeridian-safe polygon(s) | |
| """ | |
| output_format = output_format.replace("-", "").strip().lower() | |
| coords_shift = copy.deepcopy(geojson["coordinates"]) | |
| shell_minx = shell_maxx = None | |
| split_meridians = set() | |
| splitter = None | |
| for ring_index, ring in enumerate(coords_shift): | |
| if len(ring) < 1: | |
| continue | |
| else: | |
| ring_minx = ring_maxx = ring[0][0] | |
| crossings = 0 | |
| for coord_index, (lon, _) in enumerate(ring[1:], start=1): | |
| lon_prev = ring[coord_index - 1][0] # [0] corresponds to longitude coordinate | |
| if check_crossing(lon, lon_prev, validate=validate): | |
| direction = math.copysign(1, lon - lon_prev) | |
| coords_shift[ring_index][coord_index][0] = lon - (direction * 360.0) | |
| crossings += 1 | |
| x_shift = coords_shift[ring_index][coord_index][0] | |
| if x_shift < ring_minx: ring_minx = x_shift | |
| if x_shift > ring_maxx: ring_maxx = x_shift | |
| # Ensure that any holes remain contained within the (translated) outer shell | |
| if (ring_index == 0): # by GeoJSON definition, first ring is the outer shell | |
| shell_minx, shell_maxx = (ring_minx, ring_maxx) | |
| elif (ring_minx < shell_minx): | |
| ring_shift = [[x + 360, y] for (x, y) in coords_shift[ring_index]] | |
| coords_shift[ring_index] = ring_shift | |
| ring_minx, ring_maxx = (x + 360 for x in (ring_minx, ring_maxx)) | |
| elif (ring_maxx > shell_maxx): | |
| ring_shift = [[x - 360, y] for (x, y) in coords_shift[ring_index]] | |
| coords_shift[ring_index] = ring_shift | |
| ring_minx, ring_maxx = (x - 360 for x in (ring_minx, ring_maxx)) | |
| if crossings: # keep track of meridians to split on | |
| if ring_minx < -180: split_meridians.add(-180) | |
| if ring_maxx > 180: split_meridians.add(180) | |
| n_splits = len(split_meridians) | |
| if n_splits > 1: | |
| raise NotImplementedError( | |
| """Splitting a Polygon by multiple meridians (MultiLineString) | |
| not supported by Shapely""" | |
| ) | |
| elif n_splits == 1: | |
| split_lon = next(iter(split_meridians)) | |
| meridian = [[split_lon, -90.0], [split_lon, 90.0]] | |
| splitter = LineString(meridian) | |
| shell, *holes = coords_shift if splitter else geojson["coordinates"] | |
| if splitter: split_polygons = split(Polygon(shell, holes), splitter) | |
| else: split_polygons = GeometryCollection([Polygon(shell, holes)]) | |
| geo_polygons = list(translate_polygons(split_polygons, output_format)) | |
| if output_format == "geometrycollection": return GeometryCollection(geo_polygons) | |
| else: return geo_polygons |
I really need this code to work for me, I have a GeoDataFrame and can easily obtain a GeoJSON of it. I took the GeoJSON of a single polygon and passed it in, that is when I got the error above.
With a GeoDataFrame you need to use gdf.geo_interface['features'][index]['geometry'] . The GeoJSON like object that GeoPandas returns by default is a feature collection which is not what this function wants. You need to pass the geometry associated with a single polygon at a time. So you'll have to loop over all the features in the dictionary, grab their geometries independently, and then pass them to this method. By the way, you might also need to cast the coordinates in the returned geometry dictionary to a list instead of a tuple.
ah sorry, i completely missed @tony-zeidan comments! (haven't checked my github in a while)
+1 to @jgodwin's answer. the functions i wrote were intended for the GeoJSON standard (using Python dictionaries and lists), which might be different to the Python objects that GeoPandas produces. however, conversions should be quite doable. feel free to share an example of the problem you're facing on Google Colab and i'll be happy to take a look as well!
Thanks @PawaritL this is very useful code. I am just learning about antimeridian crossings myself. It's impressive how you dealt with holes etc.
I captured the above gists and started writing some python packaging, and some pytests. Also I created 1 other issue I noticed.
https://github.com/guidorice/antimeridian_splitter
Pawarit, if you want to take ownership of the repository, or add a License file, please let me know, glad to hand it to you.
Hi @guidorice, thanks so much for catching that!
It was a typo when I tried to make a quick edit directly on the gist. sorry about that
I've fixed it now and would be happy to discuss the issue you saw on geojson.io
AttributeError: 'str' object has no attribute 'bounds'is an error I got while usingsplit_polygons().