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's_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 <jswhit@...146...>
To: KURT PETERS <peterskurt@...1954...>
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 -0600KURT 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,
KurtKurt: 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 : Jeffrey.S.Whitaker@...259...
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 : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg