[Matplotlib-users] query abuot plotting polygons using a basemap projection

[Gary R. <gr...@bi...>]

I'm plotting a coverage map of a sphere using the Mollweide plot in
basemap. The attachment is an example that is produced by sending an
array of polygons (one polygon per row described as four corners, one
per column) described using polar (theta) and azimuthal (phi) angles to
the following function. As a kludge, I discard any polygons that cross
the map boundary, but this produces artefacts and it would be better to
subdivide these and keep the parts. I was wondering whether there's a
function I missed that allows me to add polygons and performs the split
across the map boundary.

Gary R.

def Mollweide(theta, phi):
      def combinations(iterable, r):
          ''' Python 2.6 itertools function'''
          # combinations('ABCD', 2) --> AB AC AD BC BD CD
          # combinations(range(4), 3) --> 012 013 023 123
          pool = tuple(iterable)
          n = len(pool)
          if r > n:
              return
          indices = range(r)
          yield tuple(pool[i] for i in indices)
          while True:
              for i in reversed(range(r)):
                  if indices[i] != i + n - r:
                      break
              else:
                  return
              indices[i] += 1
              for j in range(i+1, r):
                  indices[j] = indices[j-1] + 1
              yield tuple(pool[i] for i in indices)


      def boundary_crossed(pts):
          crossed = False
          for c in combinations(pts, 2):
              if abs(c[0]-c[1])>180:
                  crossed = True
                  break
          return crossed

      # Make Mollweide plot
      m = Basemap(projection='moll', lon_0=0, resolution='c')

      # draw the edge of the map projection region (the projection limb)
      m.drawmapboundary()
      # draw lat/lon grid lines every 30 degrees.
      m.drawmeridians(np.arange(0,360,30), dashes=[10,0])
      m.drawparallels(np.arange(-90,90,30), dashes=[10,0])

      ax = plt.gca()                          # get current axes instance
      for i in range(theta.shape[0]):
          pts = np.vstack((theta[i], phi[i])).T
          if boundary_crossed(pts[:,1]):
              continue        # skip polys that cross the map boundary

          polypts = [m(pt[1], pt[0]) for pt in pts]
          poly = Polygon(polypts, facecolor="b", edgecolor="None",
                         alpha=0.5)
          ax.add_patch(poly)