Re: [Matplotlib-users] scale a circle properly (not from shapefile)

[KURT P. <pet...@ms...>]

Thanks, that's exactly what I would like to do. I'll take a look.
Regards,
Kurt

----Original Message Follows----
From: Jeff Whitaker <js...@fa...>
To: KURT PETERS <pet...@ms...>
CC: mat...@li...
Subject: Re: [Matplotlib-users] scale a circle properly (not from shapefile)
Date: Fri, 11 Jul 2008 06:22:08 -0600

KURT PETERS wrote:
>I am trying to do something similar to the plot_tissot.py example, but am 
>having some problems.
>
>   I would like to project a group of circles onto a map projection.  Below 
>is the code I developed, which doesn't work because I get the error:
>==========ERROR ====
>   File "C:\Python25\Lib\site-packages\matplotlib\path.py", line 127, in 
>__init__
>     assert vertices.ndim == 2
>AssertionError
>==========
>
>====CODE =========
>m = Basemap(llcrnrlon=-180,llcrnrlat=-80,urcrnrlon=180,urcrnrlat=80, 
>projection='cyl')
>shp_info = m.readshapefile(r'C:\Documents and Settings\kpeters\My 
>Documents\basemap-0.99\examples\tissot','tissot',drawbounds=True)
>ax = plt.gca()
>coords = 
>[(116,45),(104,41),(98,37),(88,30),(78,25),(116,-45),(104,-41),(98,-37),(88,-30),(78,-25)]
>
>for lon1,lat1 in coords:
>     newverts   = []
>     circle = Circle((lon1,lat1),radius=10, facecolor='green')
>#    trans = circle.get_patch_transform()
>     path = circle.get_path()
>     #for jj in path.iter_segments():   #looks like the iterator is 
>broken???
>     for jj in path.vertices:
>         verts1, verts2 = jj;
>         newverts.append(m(verts1,verts2))
>     print newverts
>     p = PolyCollection(newverts, facecolor='green', zorder = 10)
>     ax.add_collection(p)
>==END CODE==
>
>Is this a logical/best way to get circles properly projected, or is there a 
>better way?
>
>I looked at "transform_vector" but I'm not too sure what the uin and vin 
>do.   Is there a transform in basemaps that could be passed to a path like 
>in this thread: "Re: [Matplotlib-users] Drawing filled circles (discs)":
>"circle = CirclePolygon((x1,y1), r, resolution)
>trans = circle.get_patch_transform()
>path = circle.get_path()
>transpath = path.transformed(trans)"
>
>It should be noted that I also tried:
>===code dif===
>for lon1,lat1 in coords:
>     newverts   = []
>     circle = Circle((lon1,lat1),radius=10, facecolor='green')
>     path = circle.get_path()
>     #for jj in path.iter_segments():   #looks like the iterator is 
>broken???
>     for jj in path.vertices:
>         verts1, verts2 = jj;
>         newverts.append(m(verts1,verts2))
>     print newverts
>#    newcircle = Circle(m(lon1,lat1),radius=10, facecolor='green')
>     p = Polygon(newverts, facecolor='green', zorder = 10)
>     ax.add_patch(p)
>===========
>but that doesn't seem to display anything (I suspect the right radius isn't 
>being used).  Note, that the "newcircle" line that is commented out, puts 
>circles on the map, they're just not transformed right.
>
>Regards,
>Kurt
>
>

Kurt:  Sounds like what you want is to create a set of points that is 
equidistant on the surface of the earth from a given central point.  I've 
cobbled something together to do this using the Geod class that is part of 
the pyproj module (http://code.google.com/p/pyproj) used by basemap.  Here 
it is:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import pyproj
from matplotlib.patches import Polygon

# Tissot's Indicatrix (http://en.wikipedia.org/wiki/Tissot's_Indicatrix).
# These diagrams illustrate the distortion inherent in all map projections.
# In conformal projections, where angles are conserved around every 
location,
# the Tissot's indicatrix are all circles, with varying sizes. In equal-area
# projections, where area proportions between objects are conserved, the
# Tissot's indicatrix have all unit area, although their shapes and
# orientations vary with location.

class Basemap2(Basemap):
    def circle(self,lon_0,lat_0,radius_deg,npts):
        # create list of npts lon,lat pairs that are equidistant on the
        # surface of the earth from central point lon_0,lat_0
        # and has radius along lon_0 of radius_deg degrees of latitude.
        # uses the WGS84 ellipsoid (a=6378137.0,b=6356752.3142).
        # points are transformed to map projection coordinates.
        g = pyproj.Geod(ellps='WGS84')
        az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
        seg = []
        delaz = 360./npts
        az = az12
        for n in range(npts):
            az = az+delaz
            lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
            seg.append(m(lon,lat))
        return seg

# create mercator Basemap instance with WGS84 ellipsoid.
m = Basemap2(llcrnrlon=-180,llcrnrlat=-70,urcrnrlon=180,urcrnrlat=70,
             projection='merc',lat_ts=20,rsphere=(6378137.0,6356752.3142))
ax = plt.gca()
# draw "circles" at specified longitudes and latitudes.
for parallel in range(-60,61,30):
    for meridian in range(-165,166,30):
        seg = m.circle(meridian,parallel,6,101)
        poly = Polygon(seg,facecolor='green',zorder=10)
        ax.add_patch(poly)
# draw meridians and parallels.
m.drawparallels(np.arange(-90,91,30),labels=[1,0,0,0])
m.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])
m.drawcoastlines()
m.fillcontinents()
plt.title('Tissot Diagram - Mercator')
plt.show()

Let us know if you have questions about what is going on here.

-Jeff



--
Jeffrey S. Whitaker         Phone : (303)497-6313
NOAA/OAR/CDC  R/PSD1        FAX   : (303)497-6449
325 Broadway                Boulder, CO, USA 80305-3328