[Matplotlib-users] Not scaling properly? Am I doing something wrong? From: KURT PETERS - 2008-03-27 14:29 ```I'm trying what I thought was a simple test and getting "bad" results. I am taking some lat long coords, and feeding it into a map. The conversion is not giving "real" values that can be plotted on a map (and actually produces an error when I use annotate). I'm including the simple code and the output: >>>>>>>>>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='ortho',lat_1=33,lat_2=45,lon_0=-95,lat_0=40,resolution='c') #COS #D + M/60 + S/3600 COSLat=38+56.0/60.0+0.013 COSLon=-1*(104+48.0/60.0) WASHLat=38+53.0/60.0+23.0/3600.0 WASHLon=-1*(77+32.0/3600.0) print COSLat x, y = m1(COSLat,COSLon) print 'x =%f, y=%f' % (x,y) m1.plot([x],[y],'ko') ax=p.gca() ax.annotate('COS1',(COSLat,COSLon)) #ax.annotate('COS2',(x,y)) ax.annotate('Wash1',(WASHLat,WASHLon)) x, y = m1(WASHLat,WASHLon) #ax.annotate('Wash2',(x,y)) 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('full resolution') p.show() <<<<<<<<>>>>>>>>>> 38.9463333333 x =1000000000000000000000000000000.000000, y=1000000000000000000000000000000.000000 Regards, Kurt ```
 Re: [Matplotlib-users] Not scaling properly? Am I doing something wrong? From: Jeff Whitaker - 2008-03-27 14:46 ```KURT PETERS wrote: > I'm trying what I thought was a simple test and getting "bad" results. I am > taking some lat long coords, and feeding it into a map. The conversion is > not giving "real" values that can be plotted on a map (and actually produces > an error when I use annotate). > I'm including the simple code and the output: > >>>>>>>>>> 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='ortho',lat_1=33,lat_2=45,lon_0=-95,lat_0=40,resolution='c') > > > #COS > #D + M/60 + S/3600 > COSLat=38+56.0/60.0+0.013 > COSLon=-1*(104+48.0/60.0) > WASHLat=38+53.0/60.0+23.0/3600.0 > WASHLon=-1*(77+32.0/3600.0) > > print COSLat > x, y = m1(COSLat,COSLon) > print 'x =%f, y=%f' % (x,y) > > m1.plot([x],[y],'ko') > ax=p.gca() > ax.annotate('COS1',(COSLat,COSLon)) > #ax.annotate('COS2',(x,y)) > ax.annotate('Wash1',(WASHLat,WASHLon)) > x, y = m1(WASHLat,WASHLon) > #ax.annotate('Wash2',(x,y)) > > 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('full resolution') > p.show() > <<<<<<<<>>>>>>>>>> > 38.9463333333 > x =1000000000000000000000000000000.000000, > y=1000000000000000000000000000000.000000 > > Regards, > Kurt > > Kurt: If you want the Lambert conformal projection, you should use projection='lcc', not 'ortho'. Nevertheless, your example works for me if I change the order of the arguments passed to the Basemap instance to x, y = m1(COSLon,COSLat) x, y = m1(WASHLon,WASHLat) Note, longitude goes first. Also, if you want 'full resolution' coastlines, use resolution='h'. -Jeff -- Jeffrey S. Whitaker Phone : (303)497-6313 NOAA/OAR/CDC R/PSD1 FAX : (303)497-6449 325 Broadway Boulder, CO, USA 80305-3328 ```
 [Matplotlib-users] Basemaps - shapefile import/display for points From: KURT PETERS - 2008-03-27 14:59 ```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 ```
 Re: [Matplotlib-users] Basemaps - shapefile import/display for points From: KURT PETERS - 2008-03-27 15:16 ```And, before someone asks, "Why are you using "h" and this line: h= [seg[0]*0.000278,seg[1]*0.000278] ax.annotate(seqnum[nshape],h)", I was using this, instead, but tried to experiment with things to try to make things work right: ax.annotate(seqnum[nshape],seg) . I usually don't "give up" and post a question unless I've tried a myriad of things, and unfortunately those things sometimes show up in the example code. Regards, Kurt ```
 Re: [Matplotlib-users] Basemaps - shapefile import/display for points From: Jeff Whitaker - 2008-03-31 12:28 ```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 ```