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 <jswhit@...146...>
To: KURT PETERS <peterskurt@...1954...>
CC: matplotlib-users@lists.sourceforge.net
Subject: Re: [Matplotlib-users] Basemaps - shapefile import/display for points
Date: Mon, 31 Mar 2008 06:27:50 -0600KURT 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#chprefWhile 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,999Here'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+=1print 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,
KurtKurt: 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 : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory