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 <[EMAIL PROTECTED]> > To: KURT PETERS <[EMAIL PROTECTED]> > CC: matplotlib-users@lists.sourceforge.net > Subject: Re: [Matplotlib-users] Basemaps - shapefile import/display > for points > Date: Mon, 31 Mar 2008 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 : [EMAIL PROTECTED] 325 Broadway Office : Skaggs Research Cntr 1D-124 Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg ------------------------------------------------------------------------- Check out the new SourceForge.net Marketplace. It's the best place to buy or sell services for just about anything Open Source. http://ad.doubleclick.net/clk;164216239;13503038;w?http://sf.net/marketplace _______________________________________________ Matplotlib-users mailing list Matplotlib-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/matplotlib-users
