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