query abuot plotting polygons using a basemap projection

[Gary Ruben-2]

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)


------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

[Jeff Whitaker]

Gary Ruben wrote:

> 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.

Gary:  You might be able to use the _geoslib module to compute the
intersections of those polygons with the map boundary.  I do a similar
thing with the coastline polygons in the _readboundarydata function.    
The _boundarypolyll and _boundarypolyxy instance variables have the
vertices of the map projection region polygons in lat/lon and projection
coords.  You could do somethig like this:

            from mpl_toolkits.basemap import _geoslib
            poly = _geoslib.Polygon(b) # a geos Polygon instance
describing your polygon)
            b = self._boundarypolyxy.boundary
            bx = b[:,0]; by= b[:,1]
            boundarypoly = _geoslib.Polygon(b) # a geos Polygon instance
describing the map region
            if poly.intersects(boundarypoly):
                    geoms = poly.intersection(boundarypoly)
                    polygons = [] # polygon intersections to plot.
                    for psub in geoms:
                            b = psub.boundary # boundary of an intersection
                            polygons.append(zip(b[:,0],b[:,1]))

-Jeff

>
> 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)
>
> ------------------------------------------------------------------------
>
> ------------------------------------------------------------------------
>
> ------------------------------------------------------------------------------
> Come build with us! The BlackBerry(R) Developer Conference in SF, CA
> is the only developer event you need to attend this year. Jumpstart your
> developing skills, take BlackBerry mobile applications to market and stay
> ahead of the curve. Join us from November 9 - 12, 2009. Register now!
> http://p.sf.net/sfu/devconference
> ------------------------------------------------------------------------
>
> _______________________________________________
> Matplotlib-users mailing list
> Matplotlib-users@...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>  

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

[Gary Ruben-2]

Thank you Jeff. I'll try out this solution.

Gary.

Jeff Whitaker wrote:
<snip>

> Gary:  You might be able to use the _geoslib module to compute the
> intersections of those polygons with the map boundary.  I do a similar
> thing with the coastline polygons in the _readboundarydata function.    
> The _boundarypolyll and _boundarypolyxy instance variables have the
> vertices of the map projection region polygons in lat/lon and projection
> coords.  You could do somethig like this:
>
>            from mpl_toolkits.basemap import _geoslib
>            poly = _geoslib.Polygon(b) # a geos Polygon instance
> describing your polygon)
>            b = self._boundarypolyxy.boundary
>            bx = b[:,0]; by= b[:,1]
>            boundarypoly = _geoslib.Polygon(b) # a geos Polygon instance
> describing the map region
>            if poly.intersects(boundarypoly):
>                    geoms = poly.intersection(boundarypoly)
>                    polygons = [] # polygon intersections to plot.
>                    for psub in geoms:
>                            b = psub.boundary # boundary of an intersection
>                            polygons.append(zip(b[:,0],b[:,1]))
>
> -Jeff

<snip>

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

[Gary Ruben-2]

Hi Jeff,

I finally had a chance to try this. I can't get it to work but I think
I'm close - for some reason, the way I'm creating the geos polygons
seems to always intersect the boundary polygon. It's hard to think of a
good minimal example for this so I've attached an example that
illustrates the problem - it tries to plot an icosahedron on the
Mollweide plot.

Gary R.

Jeff Whitaker wrote:

> Gary Ruben wrote:
>> 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.
>
> Gary:  You might be able to use the _geoslib module to compute the
> intersections of those polygons with the map boundary.  I do a similar
> thing with the coastline polygons in the _readboundarydata function.    
> The _boundarypolyll and _boundarypolyxy instance variables have the
> vertices of the map projection region polygons in lat/lon and projection
> coords.  You could do somethig like this:
>
>            from mpl_toolkits.basemap import _geoslib
>            poly = _geoslib.Polygon(b) # a geos Polygon instance
> describing your polygon)
>            b = self._boundarypolyxy.boundary
>            bx = b[:,0]; by= b[:,1]
>            boundarypoly = _geoslib.Polygon(b) # a geos Polygon instance
> describing the map region
>            if poly.intersects(boundarypoly):
>                    geoms = poly.intersection(boundarypoly)
>                    polygons = [] # polygon intersections to plot.
>                    for psub in geoms:
>                            b = psub.boundary # boundary of an intersection
>                            polygons.append(zip(b[:,0],b[:,1]))
>
> -Jeff

from mpl_toolkits.basemap import Basemap, _geoslib
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import numpy as np
from numpy import pi

icosahedron = \
[[0.53,0.,-0.53,0.53,-0.53,0.,0.53,-0.53,0.,0.53,0.,
  -0.53,0.85,0.53,0.85,0.85,0.85,0.53,-0.85,-0.53,-0.85,
  -0.85,-0.85,-0.53,0.,0.85,0.,0.,0.,-0.85,0.,0.,0.85,
  0.,-0.85,0.,0.,0.53,0.85,0.,0.85,0.53,0.53,0.,0.85,
  0.53,0.85,0.,-0.53,0.,-0.85,-0.53,-0.85,0.,-0.53,
  -0.85,0.,-0.53,0.,-0.85],
 [0.,0.85,0.,0.,0.,-0.85,0.,0.,0.85,0.,-0.85,0.,0.53,
  0.,-0.53,0.53,-0.53,0.,-0.53,0.,0.53,-0.53,0.53,0.,
  0.85,0.53,0.85,0.85,0.85,0.53,-0.85,-0.85,-0.53,-0.85,
  -0.53,-0.85,0.85,0.,0.53,0.85,0.53,0.,0.,-0.85,-0.53,
  0.,-0.53,-0.85,0.,0.85,0.53,0.,0.53,0.85,0.,-0.53,
  -0.85,0.,-0.85,-0.53],
 [0.85,0.53,0.85,0.85,0.85,0.53,-0.85,-0.85,-0.53,-0.85,
  -0.53,-0.85,0.,0.85,0.,0.,0.,-0.85,0.,0.85,0.,0.,0.,
  -0.85,0.53,0.,-0.53,0.53,-0.53,0.,0.53,-0.53,0.,0.53,
  0.,-0.53,0.53,0.85,0.,-0.53,0.,-0.85,0.85,0.53,0.,
  -0.85,0.,-0.53,0.85,0.53,0.,-0.85,0.,-0.53,0.85,0.,
  0.53,-0.85,-0.53,0.]]

icosahedron1 = \
[[0.53, 0.  ,-0.53, 0.53,-0.53, 0.  ],
 [0.  , 0.85, 0.   ,0.  , 0.  ,-0.85],
 [0.85, 0.53, 0.85 ,0.85, 0.85, 0.53]]


def pointsOnSphere():
    x,y,z = np.array(icosahedron)/0.851
    return 180/pi*np.arcsin(z), 180/pi*np.arctan2(y,x)


if __name__=='__main__':
    if 0:
        from enthought.mayavi import mlab
        x,y,z = icosahedron
        sphere = mlab.triangular_mesh(x, y, z, \
            np.arange(len(x)).reshape(-1,3), representation = 'wireframe')
        mlab.show()
        raise SystemExit
           
    # 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()

    theta, phi = pointsOnSphere()
    theta.shape = phi.shape = (-1,3)
    print theta.shape[0], 'polys'

    ax = plt.gca()                      # get current axes instance
    # create a geos Polygon instance describing the map region
    boundarypoly = _geoslib.Polygon(m._boundarypolyxy.boundary)
    for i in range(theta.shape[0]):
        pts = np.vstack((phi[i], theta[i])).T
        polypts = np.array([m(*pt) for pt in pts])  # to projection coords
        poly = _geoslib.Polygon(polypts)             # geos polygon for testing
        if poly.intersects(boundarypoly):
            for psub in poly.intersection(boundarypoly):
                b = psub.boundary           # boundary of an intersection
                polypts = zip(b[:,0],b[:,1])
                c = (1,) + tuple(np.random.random(2))   # warm colour
                poly = Polygon(polypts, facecolor=c, edgecolor=c)
                ax.add_patch(poly)
        else:
            c = tuple(np.random.random(2)) + (1,)       # cool colour
            poly = Polygon(polypts, facecolor=c, edgecolor=c)
            ax.add_patch(poly)

    plt.show()
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now.  http://p.sf.net/sfu/bobj-july
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users