Ordering polygon coordinates for plotting

Question

I have a model grid composed of many cells for which I would like to plot a shaded polygon on a matplotlib basemap.

Using pyproj, I first projected the points, before creating a polygon using shapely.geometry's Polygon class to extract the grid's exterior coordinates from. I then revert them back to WGS84 for passing to my plotting function:

grid_x_mesh, grid_y_mesh = pyproj.transform(wgs84, nplaea, grid_lons, grid_lats)
grid_x = grid_x_mesh.ravel()
grid_y = grid_y_mesh.ravel()
grid_poly = Polygon(zip(grid_x, grid_y))
grid_x, grid_y = grid_poly.exterior.coords.xy
grid_plons, grid_plats = pyproj.transform(nplaea, wgs84, grid_x, grid_y)

Then, using the matplotlib.basemap method, I projected the WSG84 coordinates to the map projection (nplaea in this case) and

grid_poly_x, grid_poly_y = m(grid_plons, grid_plats)
grid_poly_xy = zip(grid_poly_x, grid_poly_y)
grid_poly = Polygon(grid_poly_xy, facecolor='red', alpha=0.4)
plt.gca().add_patch(grid_poly)

When attempting to do so, I am getting a criss-cross pattern, which I assume has to do the ordering of the coordinates that I supplied to the polygon function.

I would think this has to do with either how I extracted the exterior coordinates, or just the ordering of the coordinate lists when I created the final polygon to plot.

Is there a clever way of ordering these properly if that is the problem?

Plotted polygon

Close-up

Answer 1

Soooo... apparently the shapely.geometry.Polygon method was drawing the polygon with all interior grid coordinates, which I realized due to the grid_plons and grid_plats having the same length as the np.ravel()'ed mesh coordinate array.

I ended up just doing a manual extraction of the external coordinates from the mesh coordinate arrays before passing them to the Polygon method (see below). Though, I imagine there may be a prettier and more general way of doing this.

Manual Extraction method:

grid_x_mesh, grid_y_mesh = pyproj.transform(wgs84, nplaea, grid_lons, grid_lats)

# The coordinates must be ordered in the order they are to be drawn
[grid_x.append(i) for i in grid_x_mesh[0,:]]
[grid_x.append(i) for i in grid_x_mesh[1:-1,-1]]
# Note that these two sides of the polygon are appended in reverse
[grid_x.append(i) for i in (grid_x_mesh[-1,:])[::-1]]
[grid_x.append(i) for i in (grid_x_mesh[1:-1,0])[::-1]]

[grid_y.append(i) for i in grid_y_mesh[0,:]]
[grid_y.append(i) for i in grid_y_mesh[1:-1,-1]]
[grid_y.append(i) for i in (grid_y_mesh[-1,:])[::-1]]
[grid_y.append(i) for i in (grid_y_mesh[1:-1,0])[::-1]]

grid_poly = Polygon(zip(grid_x, grid_y))

Answer 2

I agree there is some misarrangement of the grid coordinates. How was grid_lons created? Possibly a cleaner way to use Pyproj with Shapely geometries is to use a relatively new function shapely.ops.transform. For example:

import pyproj
from shapely.geometry import Polygon, Point
from shapely.ops import transform
from functools import partial

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'),  # WGS84 geographic
    pyproj.Proj(init='epsg:3575'))  # North Pole LAEA Europe

# Example grid cell, in long/lat
poly_g = Polygon(((5, 52), (5, 60), (15, 60), (15, 52), (5, 52)))

# Transform to projected system
poly_p = transform(project, poly_g)

The sanity of the coordinates should be preserved through the transformation (assuming that they were sane to begin with).