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 ------------------------------------------------------------------------- Sponsored by: SourceForge.net Community Choice Awards: VOTE NOW! Studies have shown that voting for your favorite open source project, along with a healthy diet, reduces your potential for chronic lameness and boredom. Vote Now at http://www.sourceforge.net/community/cca08 _______________________________________________ Matplotlib-users mailing list Matplotlib-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-users
