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