Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

orient function should work on folded polygons #1705

Closed
PBrockmann opened this issue Jan 3, 2023 · 18 comments
Closed

orient function should work on folded polygons #1705

PBrockmann opened this issue Jan 3, 2023 · 18 comments

Comments

@PBrockmann
Copy link

PBrockmann commented Jan 3, 2023

I have encountered a problem with folded polygons that despite the use of ``shapely.geometry.polygon.orient``` are not describe fully oriented counter-clockwise.

I have described the problem from the cartopy project but this issue is in fact related to shapely.

See: SciTools/cartopy#2111
Initialy from: https://stackoverflow.com/questions/74710358/how-to-detect-inner-polygons-from-a-multipolygon-shapely-object

So my question is:
How to force the orientation of a folded polygon to be oriented counter-clockwise ?

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np

fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(105)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

ax.set_global()
ax.gridlines()

transform = ccrs.Geodetic()

poly0 =  np.array([
    [-70,  10],
    [-60,  10],
    [-60,  20],
    [-70,  20],
    [-81,  20],
    [-80,  30],
    [-70,  30],
    [-70,  20],
    [-70,  10]
])
poly1 =  np.array([
    [-70,  10],
    [-60,  10],
    [-60,  20],
    [-70,  20],
    [-70,  30],
    [-80,  30],
    [-81,  20], 
    [-70,  20],
    [-70,  10]
])
polyCode = np.array([ 1, 2, 2, 2, 2, 2, 2, 2, 79])

# works with poly1 (filled correclty)
path =  matplotlib.path.Path(poly0, polyCode)

patch = matplotlib.patches.PathPatch(path, facecolor='red', edgecolor='black', transform=transform)
ax.add_patch(patch)

plt.show()

image

@jorisvandenbossche
Copy link
Member

Assuming that your poly0 and poly1 coordinates are to form Polygons as such:

import shapely
poly0 = shapely.Polygon(poly0)
poly1 = shapely.Polygon(poly1)

then the problem here is that those are not "valid" Polygons:

>>> poly0.is_valid
False
>>> poly1.is_valid
False
>>> shapely.is_valid_reason(poly0)
'Ring Self-intersection[-70 20]'
>>> shapely.is_valid_reason(poly1)
'Ring Self-intersection[-70 20]'

And many of the shapely functions are not guaranteed to work fully correctly with non-valid geometries.

If you want to represent such a polygon, you actually need a MultiPolygon:

>>> shapely.make_valid(poly0).wkt
'MULTIPOLYGON (((-60 10, -70 10, -70 20, -60 20, -60 10)), ((-80 30, -70 30, -70 20, -81 20, -80 30)))'

@idanmiara
Copy link
Contributor

idanmiara commented Jan 3, 2023

Try to embed this code in your example, it gives the same output with poly0 and poly1.

from shapely.geometry.polygon import orient, Polygon
poly_coords_oriented = orient(Polygon(poly0)).exterior.coords
path =  matplotlib.path.Path(poly_coords_oriented, polyCode)

image

I don't know cartopy, so not sure what did you expect.

@PBrockmann
Copy link
Author

PBrockmann commented Jan 3, 2023

Thanks for your answer.
Indeed I can correct the polygon and get a correct fill.

I have tried this already in my real case. And it is not working.

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random
import pickle

#-----------------------------------------
! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle

with open('./polys.pickle', "rb") as poly_file:
    polygons = pickle.load(poly_file)

#-----------------------------------------
def polygon2patch(poly, **kwargs):
    poly = shapely.geometry.polygon.orient(poly)
    path = Path.make_compound_path(
           Path(poly.exterior.coords, closed=True),
           *[Path(ring.coords) for ring in poly.interiors])
    patch = PathPatch(path, **kwargs)
    return patch

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(100)
#map_proj = ccrs.Orthographic(0, -60)
#map_proj._threshold /= 100.  # the default values is bad, users need to set them manually
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

holesNumber = []
for n,polygonA in enumerate(polygons.geoms):
    holes = []
    for m,polygonB in enumerate(polygons.geoms):
        if (m == n): continue     
        if polygonA.contains(polygonB):
            holes.append(polygonB.exterior.coords)
            holesNumber.append(m)
    if n in holesNumber: continue  # n is a hole
    polygonNew = geometry.Polygon(polygonA.exterior.coords, holes=holes)
    polygonNew = shapely.geometry.polygon.orient(polygonNew)   # Orient to be oriented counter-clockwise
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = polygon2patch(polygonNew, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)
    #ax.fill(*polygonNew.exterior.xy, transform=transform, color=random_color, lw=0.5, edgecolor="black")
    print(n, shapely.is_valid_reason(polygonNew))

ax.set_global()
ax.gridlines()
plt.show()

A lot of non valid polygons, I don't know how to correct them. in fact only 1, 2 and 11 cause trouble when filled.
If I used Path(poly.exterior.coords, closed=False), the display is correct except a bad filling but not overlaid.

@idanmiara
Copy link
Contributor

idanmiara commented Jan 4, 2023

Hi @PBrockmann

I'm not sure what's the correct output, but as @jorisvandenbossche pointed out - some of your polygons are not valid and you could try to make them valid before calling orient. you won't necessarily get a Polygon, as the output may be some collection.

You might want to try to take some ideas from this code:

import pathlib
from typing import Sequence

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry.polygon import orient

import shapely
import cartopy.crs as ccrs
import random
import pickle

#-----------------------------------------
import wget
url = 'https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle'
filename = pathlib.Path(wget.filename_from_url(url))
if not filename.is_file():
    filename = wget.download(url)
with open(filename, "rb") as poly_file:
    polygons = pickle.load(poly_file)

#-----------------------------------------
def get_paths(geometry, **kwargs) -> Sequence[Path]:
    geometry = orient(geometry)
    if geometry.geom_type in ["MultiPolygon", "MultiLineString", "GeometryCollection"]:
        lists = [get_paths(g, **kwargs) for g in geometry.geoms]
        combined = [item for sublist in lists for item in sublist]
        return combined
    elif geometry.geom_type in ["LineString"]:
        return [Path(geometry.coords)]
    elif geometry.geom_type in ["LinearRing"]:
        return [Path(geometry.coords, closed=True)]
    elif geometry.geom_type in ["Polygon"]:
        return [Path(geometry.exterior.coords, closed=True),
               *[Path(ring.coords, closed=True) for ring in geometry.interiors]]
    else:
        raise Exception(geometry)

fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(100)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

holesNumber = []
for n,polygonA in enumerate(polygons.geoms):
    polygonA = shapely.make_valid(polygonA)  # Orient to be oriented counter-clockwise
    geometry = shapely.make_valid(polygonA)
    paths = get_paths(geometry)
    compound_path = Path.make_compound_path(*paths)
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = PathPatch(compound_path, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)
    print(n, shapely.is_valid_reason(geometry))

ax.set_global()
ax.gridlines()
plt.show()

@PBrockmann
Copy link
Author

PBrockmann commented Jan 4, 2023

Thank you for your proposition but running the code I get a AttributeError: 'GeometryCollection' object has no attribute 'exterior'

The correct output is as described in the StackOverflow question.

image

With an Antartica correctly filled.

@PBrockmann
Copy link
Author

Back to a simpler test that does not work (filling incorrect) event with shapely.geometry.polygon.orient

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np
from shapely.geometry.polygon import orient, Polygon

fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(106)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

ax.set_global()
ax.gridlines()

transform = ccrs.Geodetic()

poly0Coords =  np.array([
       [-75.00049485,  18.66214045],
       [-73.00047845,  18.66214045],
       [-71.00051681,  18.66214045],
       [-69.00055717,  18.66214045],
       [-67.00053905,  18.66214045],
       [-66.99999475,  20.54502223],
       [-69.00000894,  20.54502223],
       [-70.9999975 ,  20.54502223],
       [-73.00001494,  20.54502223],
       [-75.00000269,  20.54502223],
       [-76.99999018,  20.54502223],
       [-79.00000694,  20.54502223],
       [-80.99999548,  20.54502223],
       [-83.00001164,  20.54502223],
       [-83.00001168,  22.40571069],
       [-80.99999533,  22.40572546],
       [-79.00000673,  22.40574024],
       [-76.99999109,  22.40575686],
       [-75.00000361,  22.40577533],
       [-75.00000269,  20.54502223],
       [-75.00049485,  18.66214045]
])

poly0Code = np.array([ 1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2, 2,  2,  2, 79])

poly0Coords_oriented = orient(Polygon(poly0Coords)).exterior.coords

#path =  matplotlib.path.Path(poly0Coords, poly0Code)
path =  matplotlib.path.Path(poly0Coords_oriented, poly0Code)

patch = matplotlib.patches.PathPatch(path, facecolor='red', edgecolor='black', transform=transform)
ax.add_patch(patch)

plt.show()

@idanmiara
Copy link
Contributor

Thank you for your proposition but running the code I get a AttributeError: 'GeometryCollection' object has no attribute 'exterior'

Oh, right, I've added support for this in #1690

@idanmiara
Copy link
Contributor

idanmiara commented Jan 4, 2023

Back to a simpler test that does not work (filling incorrect) event with shapely.geometry.polygon.orient

Please apply #1690 and try this:
(not sure what's the purpose of poly0Code, I've ignored it here)

from typing import Sequence

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np
from matplotlib.path import Path

import shapely
from shapely.geometry.polygon import orient, Polygon

fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(106)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

ax.set_global()
ax.gridlines()

transform = ccrs.Geodetic()

poly0Coords =  np.array([
       [-75.00049485,  18.66214045],
       [-73.00047845,  18.66214045],
       [-71.00051681,  18.66214045],
       [-69.00055717,  18.66214045],
       [-67.00053905,  18.66214045],
       [-66.99999475,  20.54502223],
       [-69.00000894,  20.54502223],
       [-70.9999975 ,  20.54502223],
       [-73.00001494,  20.54502223],
       [-75.00000269,  20.54502223],
       [-76.99999018,  20.54502223],
       [-79.00000694,  20.54502223],
       [-80.99999548,  20.54502223],
       [-83.00001164,  20.54502223],
       [-83.00001168,  22.40571069],
       [-80.99999533,  22.40572546],
       [-79.00000673,  22.40574024],
       [-76.99999109,  22.40575686],
       [-75.00000361,  22.40577533],
       [-75.00000269,  20.54502223],
       [-75.00049485,  18.66214045]
])

poly0Code = np.array([ 1,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2, 2,  2,  2, 79])

def get_paths(geometry, **kwargs) -> Sequence[Path]:
    geometry = orient(geometry)
    if geometry.geom_type in ["MultiPolygon", "MultiLineString", "GeometryCollection"]:
        lists = [get_paths(g, **kwargs) for g in geometry.geoms]
        combined = [item for sublist in lists for item in sublist]
        return combined
    elif geometry.geom_type in ["LineString"]:
        return [Path(geometry.coords)]
    elif geometry.geom_type in ["LinearRing"]:
        return [Path(geometry.coords, closed=True)]
    elif geometry.geom_type in ["Polygon"]:
        return [Path(geometry.exterior.coords, closed=True),
               *[Path(ring.coords, closed=True) for ring in geometry.interiors]]
    else:
        raise Exception(geometry)

geometry = shapely.make_valid(Polygon(poly0Coords))
paths = get_paths(geometry)
path = Path.make_compound_path(*paths)
patch = matplotlib.patches.PathPatch(path, facecolor='red', edgecolor='black', transform=transform)
ax.add_patch(patch)

plt.show()

@PBrockmann
Copy link
Author

PBrockmann commented Jan 4, 2023

Hummm no the poly0Code should be used to close properly the polygon (code 79) otherwise the filling is not correct.

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry import Polygon
import cartopy.crs as ccrs

#-----------------------------------------
def polygon2patch(poly, **kwargs):
    path = Path.make_compound_path(
           Path(poly.exterior.coords, closed=True),
           *[Path(ring.coords) for ring in poly.interiors])
    patch = PathPatch(path, **kwargs)
    return patch

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Orthographic(-80, 60)
map_proj._threshold /= 100.  # the default values is bad, users need to set them manually
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

p1 =  Polygon([(20, 20), (80, 20), (100,70), (10, 70)])
p2 =  Polygon([(40, 25), (40, 30), (45, 30), (45, 25)])
p3 =  Polygon([(50, 35), (50, 40), (55, 40)])
polygon = Polygon(p1.exterior.coords, holes=[p2.exterior.coords, p3.exterior.coords])

patch = polygon2patch(polygon, transform=transform, facecolor="red", lw=0.5, edgecolor="black")
ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

print(patch.get_path())

If Path(poly.exterior.coords, closed=True) the fill is correct. Notice the use of 79 in polygon codes.
image

If Path(poly.exterior.coords, closed=True) the fill is not correct. Notive the use of 2 in polygon codes.
image

Conclusion I need to use Path(poly.exterior.coords, closed=True) but also need to handle the case where the multi polygon is not correctly oriented.

Pffff. Quite hard problem.
Sorry I do not understand your solution.

@PBrockmann
Copy link
Author

Hi,

Tested again. Really confused in what to do.

Your proposed code does not give a correct plot.

import pathlib
from typing import Sequence

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry.polygon import orient

import shapely
import cartopy.crs as ccrs
import random
import pickle

#-----------------------------------------
import wget
url = 'https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle'
filename = pathlib.Path(wget.filename_from_url(url))
if not filename.is_file():
    filename = wget.download(url)
with open(filename, "rb") as poly_file:
    polygons = pickle.load(poly_file)

#-----------------------------------------
def get_paths(geometry, **kwargs) -> Sequence[Path]:
    if geometry.geom_type in ["MultiPolygon", "MultiLineString", "GeometryCollection"]:
        lists = [get_paths(g, **kwargs) for g in geometry.geoms]
        combined = [item for sublist in lists for item in sublist]
        return combined
    elif geometry.geom_type in ["LineString"]:
        return [Path(geometry.coords)]
    elif geometry.geom_type in ["LinearRing"]:
        return [Path(geometry.coords, closed=True)]
    elif geometry.geom_type in ["Polygon"]:
        geometry = orient(geometry)
        return [Path(geometry.exterior.coords, closed=True),
               *[Path(ring.coords, closed=True) for ring in geometry.interiors]]
    else:
        raise Exception(geometry)

fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(100)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

holesNumber = []
for n,polygonA in enumerate(polygons.geoms):
    geometry = shapely.make_valid(polygonA)
    paths = get_paths(geometry)
    compound_path = Path.make_compound_path(*paths)
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = PathPatch(compound_path, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

image

No the nearest good solution is:

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random
import pickle

#-----------------------------------------
! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys.pickle

with open('./polys.pickle', "rb") as poly_file:
    polygons = pickle.load(poly_file)

#-----------------------------------------
def polygon2patch(poly, **kwargs):
    path = Path.make_compound_path(
           Path(poly.exterior.coords),
           *[Path(ring.coords) for ring in poly.interiors])
    patch = PathPatch(path, **kwargs)
    return patch

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(-10)
#map_proj = ccrs.Orthographic(-10, -60)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

holesNumber = []
for n,polygonA in enumerate(polygons.geoms):
    holes = []
    for m,polygonB in enumerate(polygons.geoms):
        if (m == n): continue     
        if polygonA.contains(polygonB):
            holes.append(polygonB.exterior.coords)
            holesNumber.append(m)
    if n in holesNumber: continue  # n is a hole
    polygonNew = geometry.Polygon(polygonA.exterior.coords, holes=holes)
    polygonNew = shapely.geometry.polygon.orient(polygonNew)   # Orient to be oriented counter-clockwise
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = polygon2patch(polygonNew, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

But polygons near borders are not filled correctly:

image

@idanmiara
Copy link
Contributor

Hi @PBrockmann,
In my proposed code I tried to make the non-valid geometries valid before plotting them (calling orient might not be enough), I understand that you want to have some special treatments for holes which I didn't take into account in my code.

I can suggest that for pinpointing the problem it might be helpful if you provided a minimal example that reproduces the issue (maybe like your first example, if that makes sense).
Sorry, I'm not familiar enough with this issue to help further.

@PBrockmann
Copy link
Author

PBrockmann commented Jan 5, 2023

Yes, many thanks for being helpful in trying to resolve this issue.

I have tried to simplify the problem. But I am confused in the use of different parameters such as closed=True for Path structure in Matplotlib.

The big picture is to be able to draw continents and all islands with a filled color. Some polygons (extracted from a VTK process) have holes. The matplotlib Patches are the right structure to handle this I think. But I have to respect orientation and that seems non trial when you have islands connected only by a vertex. Making them valid from a shapely point of view leads to other problems (why so many different types of geometry once valid ?).

So yes confusing for a request quite simple in fact.
Thank you anyway.

@PBrockmann
Copy link
Author

PBrockmann commented Jan 5, 2023

Hi,

I have a small case extracted from my real polygons that explains the problem and that I cannot solve.
I fact the polygon is folded and orient it has no effect.

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random

#-----------------------------------------
def polygon2patch(poly, **kwargs):
    path = Path.make_compound_path(
           Path(poly.exterior.coords, closed=True),
           *[Path(ring.coords) for ring in poly.interiors])
    patch = PathPatch(path, **kwargs)
    return patch

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

#map_proj = ccrs.Robinson(-10)
map_proj = ccrs.Orthographic(-10, 60)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

polygonNew = Polygon([
       [146.38076034,  74.032728  ],
       [144.69941506,  74.96788777],
       [142.81028083,  75.88109748],
       [139.61828743,  75.42864436],
       [137.39423796,  76.29102758],
       [134.30816041,  75.7765326 ],
       [136.59781355,  74.94057655],
       [138.65328914,  74.0755147 ],
       [141.60305404,  74.53908526],
       [144.69941506,  74.96788777],
       [147.94064981,  75.35911624],
       [149.50617426,  74.4036397 ],
       [146.38076034,  74.032728  ]]
)

random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
patch = polygon2patch(polygonNew, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

image

Inpecting the coordinates show the problem:

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1, 1, 1)

polygonNew = Polygon([
       [146.38076034,  74.032728  ],
       [144.69941506,  74.96788777],
       [142.81028083,  75.88109748],
       [139.61828743,  75.42864436],
       [137.39423796,  76.29102758],
       [134.30816041,  75.7765326 ],
       [136.59781355,  74.94057655],
       [138.65328914,  74.0755147 ],
       [141.60305404,  74.53908526],
       [144.69941506,  74.96788777],
       [147.94064981,  75.35911624],
       [149.50617426,  74.4036397 ],
       [146.38076034,  74.032728  ]]
)

coords = np.array(polygonNew.exterior.coords)
ax.plot(coords[:,0], coords[:,1], marker='.')
ax.grid()
plt.show()

a

The polygon coordinates is folded. When it is fully describes anti-clockwise then the filling is correct.

So how to do this ? Correct all polygons.

@idanmiara
Copy link
Contributor

The polygon coordinates is folded. When it is fully describes anti-clockwise then the filling is correct.

So how to do this ? Correct all polygons.

If many of your polygons are not valid maybe it explains all the filling problems, so maybe run make_valid on all and see if you could figure out how to progress from there.
Note that the output of make_valid is not always a polygon (in your example it is a MultiPolygon).
I did try to run makr_valid and then orient on your inputs
But then I didn't go forward with the holes processing and the poly0code - so maybe you want to continue from there.

@PBrockmann
Copy link
Author

PBrockmann commented Mar 22, 2023

Hi,

I come back on this issue. Because folded polygons cannot be drawn on a map.

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapely
from shapely import geometry
import cartopy.crs as ccrs
import random

#-----------------------------------------
def polygon2patch(poly, **kwargs):
    path = Path.make_compound_path(
           Path(poly.exterior.coords, closed=True),                  
           *[Path(ring.coords) for ring in poly.interiors])
    patch = PathPatch(path, **kwargs)
    return patch

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(100)
#map_proj = ccrs.Orthographic(-10, 0)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

polygon1 = shapely.Polygon([(0,10),(10,10),(20,10),(20,0),(10,0),(10,10),(10,20),(0,20)])
print(polygon1.exterior.is_ccw)

polygon2 = shapely.Polygon([(0,10),(10,10),(10,0),(20,0),(20,10),(10,10),(10,20),(0,20)])
print(polygon2.exterior.is_ccw)

random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
patch = polygon2patch(polygon1, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

The first polygon is not drawn, the second is drawn.

image

The first one is folded and detected as oriented counter clock-wise. That is partially true.

How can I repair this orientation ?

@PBrockmann
Copy link
Author

PBrockmann commented Apr 13, 2023

Any news on this problem ?
I haven't been able to find a solution.

image

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import matplotlib

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(100)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

polygon1 = np.array([(0,10),(10,10),(20,10),(20,0),(10,0),(10,10),(10,20),(0,20),(0,10)])
polygon2 = np.array([(0,10),(10,10),(10,0),(20,0),(20,10),(10,10),(10,20),(0,20),(0,10)])

path =  matplotlib.path.Path(polygon1, closed=True)
patch = matplotlib.patches.PathPatch(path, facecolor='red', edgecolor='black', transform=transform)
ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

polygon1 is not drawn.
polygon2 is drawn correctly because order of its vertex follows ccw.

So how to correct polygon1 ?

@PBrockmann
Copy link
Author

PBrockmann commented Apr 20, 2023

Here is a solution I have found.

The strategy is to separate the folded polygon into multiple polygons, orient them correclty with shapely then build a patch matplotlib object.

Here is the code to correct polygon1 that is a folded polygon.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import matplotlib
import shapely

#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(100)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

# This polygon is folded --> filling is not correct and orientation cannot be done correctly
polygon1 = np.array([(0,10),(10,10),(20,10),(20,0),(10,0),(10,10),(10,20),(0,20),(0,10)])

# Generate ids of points
ids = [polygon1.tolist().index(value) for value in polygon1.tolist()]

# Find duplicate points
unique, count = np.unique(ids, return_counts=True)
dup = unique[count > 1]
print("---- Duplicate points: ", dup)

# Extract points between duplicate points
polys = []
for id in dup[::-1]:
    select = np.where(ids == id)[0]
    polys.append(ids[select[0]:select[1]+1])
    idsIndices = list(range(select[0],select[1]))
    ids = list(np.delete(ids, idsIndices))
print("Number of polygons found: ", len(polys))

# Create shapely polygons for counter clockwise orientation
for poly in polys:
    pointsPoly = []
    for id in poly:
        pointsPoly.append(polygon1[id])
        
    # Create shapely polygon and orient it properly
    polygon = shapely.geometry.Polygon(pointsPoly)
    polygon = shapely.geometry.polygon.orient(polygon)

    path =  matplotlib.path.Path(polygon.exterior.coords, closed=True)
    patch = matplotlib.patches.PathPatch(path, facecolor='red', edgecolor='black', transform=transform)
    ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

@PBrockmann
Copy link
Author

PBrockmann commented Apr 20, 2023

And to close this issue.
I also give my solution applied to the real problem I had to solve.

image

So same strategy. Separate the different self intersected polygons from their duplicate points into multiple polygons.
Then build a path with eventually holes, with option closed=True. Then a patch object.
Then no overlay, the filling of each polygons is correct whatever the projection.

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib
import shapely
import cartopy.crs as ccrs
import random
import pickle

#-----------------------------------------
! wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/polys_03.pickle

with open('./polys_03.pickle', "rb") as poly_file:
    polygons = pickle.load(poly_file)
    
#-----------------------------------------
def polygon2path(poly):
    path = Path.make_compound_path(
               Path(poly.exterior.coords, closed=True),
               *[Path(ring.coords) for ring in poly.interiors])
    return path
#-----------------------------------------
fig = plt.figure(figsize=(10,5))

map_proj = ccrs.Robinson(10)
#map_proj = ccrs.Orthographic(150, -40) 
#map_proj = ccrs.Orthographic(-10, 60)
ax = fig.add_subplot(1, 1, 1, projection=map_proj)

transform = ccrs.Geodetic()

paths = []
holesNumber = []
for n,polygonA in enumerate(polygons):
    holes = []
    for m,polygonB in enumerate(polygons):
        if (m == n): continue
        if polygonA.contains(polygonB):
            holes.append(polygonB.exterior.coords)
            holesNumber.append(m)
    if n in holesNumber: continue  # n is a hole
    polygonNew = shapely.geometry.Polygon(polygonA.exterior.coords, holes=holes)
    polygonNew = shapely.geometry.polygon.orient(polygonNew)   # Orient to be oriented counter-clockwise
    
    path = polygon2path(polygonNew)
    paths.append(path)

    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    patch = PathPatch(path, transform=transform, facecolor=random_color, lw=0.5, edgecolor="black")
    ax.add_patch(patch)

ax.set_global()
ax.gridlines()
plt.show()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants