KURT PETERS wrote:
> Thanks, that's exactly what I would like to do. I'll take a look.
> Regards,
> Kurt
>
Kurt: I just added a "tissot" method to Basemap that does this (so you
won't have to extend the Basemap class in the next version). The
plot_tissot.py example has been updated to do this. Here's the method:
def tissot(self,lon_0,lat_0,radius_deg,npts):
"""
create list of ``npts`` x,y pairs that are equidistant on the
surface of the earth from central point ``lon_0,lat_0`` and form
an ellipse with radius of ``radius_deg`` degrees of latitude along
longitude ``lon_0``.
The ellipse represents a Tissot's indicatrix
(http://en.wikipedia.org/wiki/Tissot%27s_Indicatrix),
which when drawn on a map shows the distortion
inherent in the map projection."""
g = pyproj.Geod(a=self.rmajor,b=self.rminor)
az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
seg = [self(lon_0,lat_0+radius_deg)]
delaz = 360./npts
az = az12
for n in range(npts):
az = az+delaz
# skip segments along equator (Geod can't handle equatorial
arcs)
if np.allclose(0.,lat_0) and (np.allclose(90.,az) or
np.allclose(270.,az)):
continue
else:
lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
x,y = self(lon,lat)
# add segment if it is in the map projection region.
if x < 1.e20 and y < 1.e20:
seg.append((x,y))
return seg
-Jeff
> ----Original Message Follows----
> From: Jeff Whitaker <[EMAIL PROTECTED]>
> To: KURT PETERS <[EMAIL PROTECTED]>
> CC: matplotlib-users@lists.sourceforge.net
> 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
>
>
>
> -------------------------------------------------------------------------
> 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
>
--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : [EMAIL PROTECTED]
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : [EMAIL PROTECTED]
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
-------------------------------------------------------------------------
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
