Thread: Re: [Matplotlib-users] Basemaps - shapefile import/display for points

KURT PETERS wrote:
> Thanks,
> I'll give that a try. I had seen the other example, but had a very 
> difficult time figuring out what this line does:
> x, y = zip(*m.cities)
Kurt:
See the docs for the zip built-in python function at 
http://docs.python.org/lib/built-in-funcs.html.
x,y = zip(*<list of tuples> is just the reverse of zip(x,y)
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/103702
>
> Frankly, I have google'ed possibilities, but "zip" is so 
> ubiquitous, that figuring out what it really does in THIS case is 
> difficult. Do you have a good explanation of why that's necessary? (I 
> saw nothing in the shapefile docs that talks to zipping files.
> I'm not sure why 'enumerate' doesn't work?
> I will give the annotate a try with the code you provided and see if 
> that works. Also, is there a particular reason why you chose '10' 
> for your zorder? Use of that parameter is not especially clear in the 
> documentaion - perhaps having a table with what other thing's zorders 
> are would help.
There's nothing magic about 10. It just has to be greater than 1 so the 
dots come out on top of the continent fill. If you leave out the call 
to fillcontinents, you don't need the zorder at all.
>
> As for suggestions about the shapefile doc/usability. I think it's 
> hard to handle such a multidimensional data format in a workable 
> sense. I'm getting the hang of it. It's hard to have visualization 
> of what the shapefile looks like. Perhaps some kind of auto schematic 
> (think Visio or graphviz) function would be neat to show how things 
> are mapped in the shapefile and something that tells you line, poly, 
> point, etc., in the blocks and how they map to a built-in pylib class?
The thing is, shapefiles are not really a format I use a lot. I tend to 
work on things that I actually use the most.
> If you had more of a wiki format to the documentation, I know I 
> would have modified the docs to make things clearer as I've been 
> muddling through. Perhaps making a tutorial. It's especially true 
> since I've been teaching myself about shapefiles as well.
That's a good idea. Making a tutorial has always been on my to-do list, 
but I never seem to get the time. It would be a great if someone would 
step up and contribute one.
-Jeff
>
> Kurt
> ----Original Message Follows----
> From: Jeff Whitaker <js...@fa...>
> To: KURT PETERS <pet...@ms...>
> CC: mat...@li...
> Subject: Re: [Matplotlib-users] Basemaps - shapefile import/display 
> for points
> Date: 2008年3月31日 06:27:50 -0600
>
> KURT PETERS wrote:
>> OK Jeff, Thanks for your help on the previous question - I had been 
>> playing with different projections and resolutions, so that's why the 
>> comments didn't match the actual settings in the procedure calls. 
>> Now for a "real" problem:
>>
>> I'm trying to plot the cities from this web site: 
>> http://nationalatlas.gov/metadata/citiesx020.faq.html
>> using that shapefile, which uses points, not polygons (it took a long 
>> time to figure out that difference from the example of fillstates.py).
>> http://nationalatlas.gov/atlasftp.html?openChapters=chpref#chpref
>>
>> While I think I'm loading everything and displaying everything 
>> correctly, the values are not plotting right, nor do they seem 
>> realistic.
>>
>> For instance the point values look like this (which really can't be 
>> right):
>>
>> Shape num Fairbanks, coords=(42082.855349492747, 5336578.2660309337)
>> Shape num Anchorage, coords=(-442294.67146861833, 5031412.4918638617)
>>
>> print shp_info - the second value shows to use points not polys:
>> (35432, 1, [-174.20294189453125, 17.711706161499023, 0.0, 0.0], 
>> [178.87460327148437, 71.290138244628906, 0.0, 0.0])
>> Dictionaries:
>> ['STATE_FIPS', 'NAME', 'POP_2000', 'FEATURE', 'COUNTY', 'STATE', 
>> 'FIPS', 'CITIESX020', 'FIPS55', 'DISPLAY', 'POP_RANGE']
>> STATE_FIPS = 02, NAME = Anchorage, POP_2000=260283, FEATURE = County 
>> Seat, COUNTY=Anchorage Borough, STATE=AK, FIPS=02020, CITIESX020 = 
>> 194, FIPS55=03000, DISPLAY=0, POP_RANGE=250,000 - 499,999
>>
>>
>>
>> Here's the code:
>> ===============
>> import pylab as p
>> import numpy
>> from matplotlib.toolkits.basemap import Basemap as Basemap
>> from matplotlib.colors import rgb2hex
>> from matplotlib.patches import Polygon
>>
>> # Lambert Conformal map of lower 48 states.
>> # create new figure
>> fig=p.figure()
>> m1 = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,\
>> projection='lcc',lat_1=33,lat_2=45,lon_0=-95,resolution='c')
>> shp_info = 
>> m1.readshapefile(r'C:\Python25\Lib\basemap-0.9.9.1\examples\citiesx020','states',drawbounds=True) 
>>
>>
>> ax=p.gca()
>>
>> #define SHPT_POINT 1 Points
>> #define SHPT_ARC 3 Arcs (Polylines, possible in parts)
>> #define SHPT_POLYGON 5 Polygons (possible in parts)
>> #define SHPT_MULTIPOINT 8 MultiPoint (related points)
>> print shp_info
>> print m1.states_info[0].keys()
>> seqnum={}
>> criteriatodisplay=[]
>> ii=0
>> for shapedict in m1.states_info:
>> if int(shapedict['POP_2000'])>100000:
>> #'STATE_FIPS', 'NAME', 'POP_2000', 'FEATURE', 'COUNTY', 'STATE', 
>> 'FIPS', 'CITIESX020', 'FIPS55', 'DISPLAY', 'POP_RANGE']
>> print 'STATE_FIPS = %s, NAME = %s, POP_2000=%s, FEATURE = %s, 
>> COUNTY=%s, STATE=%s, FIPS=%s, CITIESX020 = %s, FIPS55=%s, DISPLAY=%s, 
>> POP_RANGE=%s' %\
>> (str(shapedict['STATE_FIPS']), str(shapedict['NAME']), 
>> str(shapedict['POP_2000']), str(shapedict['FEATURE']), 
>> str(shapedict['COUNTY']), str(shapedict['STATE']), 
>> str(shapedict['FIPS']), str(shapedict['CITIESX020']), 
>> str(shapedict['FIPS55']), str(shapedict['DISPLAY']), 
>> str(shapedict['POP_RANGE']))
>> seqnum[shapedict['CITIESX020']]=shapedict['NAME']
>> criteriatodisplay.append(shapedict['CITIESX020'])
>> ii+=1
>>
>> print ii
>>
>> for nshape,seg in enumerate(m1.states):
>> if nshape in criteriatodisplay:
>> print 'Shape num %s, coords=%s' % (seqnum[nshape], seg)
>> h= [seg[0]*0.000278,seg[1]*0.000278]
>>
>> ax.annotate(seqnum[nshape],h)
>> m1.drawcoastlines()
>> m1.fillcontinents()
>> m1.drawcountries()
>> m1.drawstates()
>> m1.drawparallels(numpy.arange(25,65,4),labels=[1,0,0,0])
>> m1.drawmeridians(numpy.arange(-120,-40,4),labels=[0,0,0,1])
>> p.title('Test Cities')
>> p.show()
>> =============
>> Regards,
>> Kurt
>>
>>
>
> Kurt: I had no trouble plotting them with this script:
>
> m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,\
> projection='lcc',lat_1=33,lat_2=45,lon_0=-95,resolution='c')
> shp_info = m.readshapefile('citiesx020','cities')
> x, y = zip(*m.cities)
> m.drawcoastlines()
> m.drawcountries()
> m.fillcontinents()
> m.scatter(x,y,2,'b',marker='o',faceted=False,zorder=10)
> p.show()
>
> This is adapted from the plotcities.py example, which was designed for 
> point files (fillstates.py was designed for polygon files). In this 
> case, m.cities is just a list of x,y coordinates. I don't know why 
> ax.annotate wasn't working for you.
>
> I know the shapefile stuff is non-intuitive and could use a lot of 
> work. Perhaps when you can offer some suggestions for the docs, or 
> for re-designing the interface.
>
> -Jeff
>
>
> -- 
> Jeffrey S. Whitaker Phone : (303)497-6313
> NOAA/OAR/CDC R/PSD1 FAX : (303)497-6449
> 325 Broadway Boulder, CO, USA 80305-3328
>
>
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

can someone explain why scatter would work but gca.annotate would not when 
plotting data on a map (see previous posts)?
 I've also tried pylib.figtext and that doesn't work either.
Regards,
Kurt

KURT PETERS wrote:
> can someone explain why scatter would work but gca.annotate would not when 
> plotting data on a map (see previous posts)?
> I've also tried pylib.figtext and that doesn't work either.
>
> Regards,
> Kurt
>
>
> 
Kurt: This works for me (with the latest svn version of basemap and 
matplotlib). Does it not work for you?
import pylab as p
from mpl_toolkits.basemap import Basemap as Basemap # use for svn
#from matplotlib.toolkits.basemap import Basemap # use for released version
m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,\
 projection='lcc',lat_1=33,lat_2=45,lon_0=-95,resolution='c')
shp_info = m.readshapefile('citiesx020','cities')
for xy in m.cities:
 p.gca().annotate('x',xy)
m.drawcoastlines()
p.show()
-Jeff
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Jeff,
 Do you think it's possible the names or CITIESX020 variable are not being 
brought in in the right order?
 I modified my code to use scatter, and, although it looks like the dots 
are in the right place, the names aren't matching?
see code:
=========
import pylab as p
import numpy
from matplotlib.toolkits.basemap import Basemap as Basemap
from matplotlib.colors import rgb2hex
from matplotlib.patches import Polygon
# Lambert Conformal map of lower 48 states.
# create new figure
#http://nationalatlas.gov/metadata/citiesx020.faq.html
#http://nationalatlas.gov/atlasftp.html?openChapters=chpref#chpref
fig=p.figure()
m1 = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,\
 projection='lcc',lat_1=33,lat_2=45,lon_0=-95,resolution='c')
shp_info = 
m1.readshapefile(r'C:\Python25\Lib\basemap-0.9.9.1\examples\citiesx020','states',drawbounds=True)
ax=p.gca()
#define SHPT_POINT 1 Points
#define SHPT_ARC 3 Arcs (Polylines, possible in parts)
#define SHPT_POLYGON 5 Polygons (possible in parts)
#define SHPT_MULTIPOINT 8 MultiPoint (related points)
print shp_info
print m1.states_info[0].keys()
seqnum={}
criteriatodisplay=[]
names={}
ii=0
for shapedict in m1.states_info:
 if int(shapedict['POP_2000'])>100000:
#'STATE_FIPS', 'NAME', 'POP_2000', 'FEATURE', 'COUNTY', 'STATE', 'FIPS', 
'CITIESX020', 'FIPS55', 'DISPLAY', 'POP_RANGE']
 print 'STATE_FIPS = %s, NAME = %s, POP_2000=%s, FEATURE = %s, 
COUNTY=%s, STATE=%s, FIPS=%s, CITIESX020 = %s, FIPS55=%s, DISPLAY=%s, 
POP_RANGE=%s' %\
 (str(shapedict['STATE_FIPS']), str(shapedict['NAME']), 
str(shapedict['POP_2000']), str(shapedict['FEATURE']), 
str(shapedict['COUNTY']), str(shapedict['STATE']), str(shapedict['FIPS']), 
str(shapedict['CITIESX020']), str(shapedict['FIPS55']), 
str(shapedict['DISPLAY']), str(shapedict['POP_RANGE']))
 # have an index of the names
 seqnum[shapedict['CITIESX020']]=shapedict['NAME']
 criteriatodisplay.append(shapedict['CITIESX020'])
 ii+=1
print ii
#x, y = zip(*m1.states)
#print m1.states[1]
#print x[1]
#print y[1]
#ii=0
x=[]
y=[]
ii=0
for nshape,seg in enumerate(m1.states):
 if nshape in criteriatodisplay:
 x.append(seg[0])
 y.append(seg[1])
 p.text(seg[0],seg[1],seqnum[nshape],fontsize=12)
 ii+=1
 #print 'Shape num %s, coords=%s' % (seqnum[nshape], seg)
 # ax.annotate(seqnum[nshape],seg)
m1.scatter(x,y,2,'b',marker='o',faceted=False,zorder=10)
#ax.annotate(s='s',xy=(int(x),int(y)))
#p.figtext(x,y,'o',weight='heavy', size = 16)
m1.drawcoastlines()
m1.fillcontinents()
m1.drawcountries()
m1.drawstates()
m1.drawparallels(numpy.arange(25,65,4),labels=[1,0,0,0])
m1.drawmeridians(numpy.arange(-120,-40,4),labels=[0,0,0,1])
p.title('Test Cities')
p.show()
========
kurt
<snip>

KURT PETERS wrote:
> Jeff,
> Do you think it's possible the names or CITIESX020 variable are not 
> being brought in in the right order?
> I modified my code to use scatter, and, although it looks like the 
> dots are in the right place, the names aren't matching?
> see code:
> =========
> import pylab as p
> import numpy
> from matplotlib.toolkits.basemap import Basemap as Basemap
> from matplotlib.colors import rgb2hex
> from matplotlib.patches import Polygon
>
> # Lambert Conformal map of lower 48 states.
> # create new figure
> #http://nationalatlas.gov/metadata/citiesx020.faq.html
> #http://nationalatlas.gov/atlasftp.html?openChapters=chpref#chpref
> fig=p.figure()
> m1 = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,\
> projection='lcc',lat_1=33,lat_2=45,lon_0=-95,resolution='c')
> shp_info = 
> m1.readshapefile(r'C:\Python25\Lib\basemap-0.9.9.1\examples\citiesx020','states',drawbounds=True) 
>
>
> ax=p.gca()
>
> #define SHPT_POINT 1 Points
> #define SHPT_ARC 3 Arcs (Polylines, possible in parts)
> #define SHPT_POLYGON 5 Polygons (possible in parts)
> #define SHPT_MULTIPOINT 8 MultiPoint (related points)
> print shp_info
> print m1.states_info[0].keys()
> seqnum={}
> criteriatodisplay=[]
> names={}
> ii=0
> for shapedict in m1.states_info:
> if int(shapedict['POP_2000'])>100000:
> #'STATE_FIPS', 'NAME', 'POP_2000', 'FEATURE', 'COUNTY', 'STATE', 
> 'FIPS', 'CITIESX020', 'FIPS55', 'DISPLAY', 'POP_RANGE']
> print 'STATE_FIPS = %s, NAME = %s, POP_2000=%s, FEATURE = %s, 
> COUNTY=%s, STATE=%s, FIPS=%s, CITIESX020 = %s, FIPS55=%s, DISPLAY=%s, 
> POP_RANGE=%s' %\
> (str(shapedict['STATE_FIPS']), str(shapedict['NAME']), 
> str(shapedict['POP_2000']), str(shapedict['FEATURE']), 
> str(shapedict['COUNTY']), str(shapedict['STATE']), 
> str(shapedict['FIPS']), str(shapedict['CITIESX020']), 
> str(shapedict['FIPS55']), str(shapedict['DISPLAY']), 
> str(shapedict['POP_RANGE']))
> # have an index of the names
> seqnum[shapedict['CITIESX020']]=shapedict['NAME']
> criteriatodisplay.append(shapedict['CITIESX020'])
>
> ii+=1
>
> print ii
> #x, y = zip(*m1.states)
> #print m1.states[1]
> #print x[1]
> #print y[1]
> #ii=0
> x=[]
> y=[]
> ii=0
> for nshape,seg in enumerate(m1.states):
> if nshape in criteriatodisplay:
> x.append(seg[0])
> y.append(seg[1])
> p.text(seg[0],seg[1],seqnum[nshape],fontsize=12)
> ii+=1
>
> #print 'Shape num %s, coords=%s' % (seqnum[nshape], seg)
> # ax.annotate(seqnum[nshape],seg)
> m1.scatter(x,y,2,'b',marker='o',faceted=False,zorder=10)
>
> #ax.annotate(s='s',xy=(int(x),int(y)))
> #p.figtext(x,y,'o',weight='heavy', size = 16)
> m1.drawcoastlines()
> m1.fillcontinents()
> m1.drawcountries()
> m1.drawstates()
> m1.drawparallels(numpy.arange(25,65,4),labels=[1,0,0,0])
> m1.drawmeridians(numpy.arange(-120,-40,4),labels=[0,0,0,1])
> p.title('Test Cities')
> p.show()
> ========
> kurt
>
> <snip>
>
>
Kurt: If I look at the least entry in m.states_info (Newport, VT), the 
corresponding x,y location in m.states corresponds to 44.93N, -72.21W, 
which looks about right to me. I guess I'm still not clear on what the 
issue is. Could you distill your example code down to something very 
simple that clearly demonstrates the problem?
-Jeff
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

OK, I figured out what I was doing wrong. Quite simple really:
I assumed the CITIESx020 value was the same as the "enumerated" value - 
which wasn't the case. Thus the wrong city name was being assigned to some 
correct city's location. Once I figured that out, all seems better.
One thing I DID notice though, which seems to be a bug in the rendering tool 
is this:
When I zoom into an area using the toolbar at the bottom, the annotations 
appear outside the "zoomed in region", ie. the rectangular area zoomed in 
to.
Regards,
Kurt
----Original Message Follows----
From: Jeff Whitaker <js...@fa...>
To: KURT PETERS <pet...@ms...>
CC: mat...@li...
Subject: Re: [Matplotlib-users] Basemaps - shapefile import/display for 
points
Date: 2008年4月01日 15:54:48 -0600
KURT PETERS wrote:
>Jeff,
> Do you think it's possible the names or CITIESX020 variable are not being 
>brought in in the right order?
> I modified my code to use scatter, and, although it looks like the dots 
>are in the right place, the names aren't matching?
>see code:
>=========
>import pylab as p
>import numpy
>from matplotlib.toolkits.basemap import Basemap as Basemap
>from matplotlib.colors import rgb2hex
>from matplotlib.patches import Polygon
>
># Lambert Conformal map of lower 48 states.
># create new figure
>#http://nationalatlas.gov/metadata/citiesx020.faq.html
>#http://nationalatlas.gov/atlasftp.html?openChapters=chpref#chpref
>fig=p.figure()
>m1 = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,\
> projection='lcc',lat_1=33,lat_2=45,lon_0=-95,resolution='c')
>shp_info = 
>m1.readshapefile(r'C:\Python25\Lib\basemap-0.9.9.1\examples\citiesx020','states',drawbounds=True)
>
>
>ax=p.gca()
>
>#define SHPT_POINT 1 Points
>#define SHPT_ARC 3 Arcs (Polylines, possible in parts)
>#define SHPT_POLYGON 5 Polygons (possible in parts)
>#define SHPT_MULTIPOINT 8 MultiPoint (related points)
>print shp_info
>print m1.states_info[0].keys()
>seqnum={}
>criteriatodisplay=[]
>names={}
>ii=0
>for shapedict in m1.states_info:
> if int(shapedict['POP_2000'])>100000:
>#'STATE_FIPS', 'NAME', 'POP_2000', 'FEATURE', 'COUNTY', 'STATE', 'FIPS', 
>'CITIESX020', 'FIPS55', 'DISPLAY', 'POP_RANGE']
> print 'STATE_FIPS = %s, NAME = %s, POP_2000=%s, FEATURE = %s, 
>COUNTY=%s, STATE=%s, FIPS=%s, CITIESX020 = %s, FIPS55=%s, DISPLAY=%s, 
>POP_RANGE=%s' %\
> (str(shapedict['STATE_FIPS']), str(shapedict['NAME']), 
>str(shapedict['POP_2000']), str(shapedict['FEATURE']), 
>str(shapedict['COUNTY']), str(shapedict['STATE']), str(shapedict['FIPS']), 
>str(shapedict['CITIESX020']), str(shapedict['FIPS55']), 
>str(shapedict['DISPLAY']), str(shapedict['POP_RANGE']))
> # have an index of the names
> seqnum[shapedict['CITIESX020']]=shapedict['NAME']
> criteriatodisplay.append(shapedict['CITIESX020'])
>
> ii+=1
>
>print ii
>#x, y = zip(*m1.states)
>#print m1.states[1]
>#print x[1]
>#print y[1]
>#ii=0
>x=[]
>y=[]
>ii=0
>for nshape,seg in enumerate(m1.states):
> if nshape in criteriatodisplay:
> x.append(seg[0])
> y.append(seg[1])
> p.text(seg[0],seg[1],seqnum[nshape],fontsize=12)
> ii+=1
>
> #print 'Shape num %s, coords=%s' % (seqnum[nshape], seg)
> # ax.annotate(seqnum[nshape],seg)
>m1.scatter(x,y,2,'b',marker='o',faceted=False,zorder=10)
>
>#ax.annotate(s='s',xy=(int(x),int(y)))
>#p.figtext(x,y,'o',weight='heavy', size = 16)
>m1.drawcoastlines()
>m1.fillcontinents()
>m1.drawcountries()
>m1.drawstates()
>m1.drawparallels(numpy.arange(25,65,4),labels=[1,0,0,0])
>m1.drawmeridians(numpy.arange(-120,-40,4),labels=[0,0,0,1])
>p.title('Test Cities')
>p.show()
>========
>kurt
>
><snip>
>
>
Kurt: If I look at the least entry in m.states_info (Newport, VT), the 
corresponding x,y location in m.states corresponds to 44.93N, -72.21W, which 
looks about right to me. I guess I'm still not clear on what the issue is. 
Could you distill your example code down to something very simple that 
clearly demonstrates the problem?
-Jeff
--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

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

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: 2008年7月11日 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