Python Basemap.readshapefile方法代码示例

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import readshapefile [as 别名]
"""
draw Atlantic Hurricane Tracks for storms that reached Cat 4 or 5.
part of the track for which storm is cat 4 or 5 is shown red.
ESRI shapefile data from http://www.nationalatlas.gov/atlasftp.html
"""
import pylab as p
from matplotlib.toolkits.basemap import Basemap as Basemap
# Lambert Conformal Conic map.
m = Basemap(llcrnrlon=-100.,llcrnrlat=0.,urcrnrlon=-20.,urcrnrlat=57.,
            projection='lcc',lat_1=20.,lat_2=40.,lon_0=-60.)
# make sure map has right aspect ratio.
fig=p.figure(figsize=(8,m.aspect*8))
fig.add_axes([0.1,0.1,0.8,0.8])
# read shapefile.
shp_info = m.readshapefile('huralll020','hurrtracks',drawbounds=False)
print shp_info
# find names of storms that reached Cat 4.
names = []
for shapedict in m.hurrtracks_info:
    cat = shapedict['CATEGORY']
    name = shapedict['NAME']
    if cat in ['H4','H5'] and name not in names:
        if name != 'NOT NAMED':  names.append(name)
print names
print len(names)
# plot tracks of those storms.
for shapedict,shape in zip(m.hurrtracks_info,m.hurrtracks):
    name = shapedict['NAME']
    cat = shapedict['CATEGORY']
    if name in names:
        xx,yy = zip(*shape)

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import readshapefile [as 别名]
from matplotlib.toolkits.basemap import Basemap as Basemap
from matplotlib.collections import LineCollection
from matplotlib.colors import rgb2hex
import random

# requires pyshapelib from Thuban (http://thuban.intevation.org/).
# cd to libraries/pyshapelib in Thuban source distribution, run
# 'python setup.py install'.

# Lambert Conformal map of lower 48 states.
m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
fig=p.figure(figsize=(8,m.aspect*8))
fig.add_axes([0.1,0.1,0.8,0.8])
# draw climate division boundaries.
shp_info = m.readshapefile('divisions','climdivs',drawbounds=True)
print shp_info
# make sure the shapefile has polygons (and not just lines).
if shp_info[1] != 5:
    print 'warning: shapefile does not contain polygons'
# choose a color for each climate division (randomly).
colors={}
divnames=[]
print m.climdivs_info[0].keys()
for shapedict in m.climdivs_info:
    divname = shapedict['ST']+repr(shapedict['DIV'])
    colors[divname] = (random.uniform(0,1),random.uniform(0,1),random.uniform(0,1))
    divnames.append(divname)
# cycle through climate divnames, color each one.
for nshape,seg in enumerate(m.climdivs):
    xx,yy = zip(*seg)

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import readshapefile [as 别名]
lon_150=[]
lat_150=[]
for xy in verts:
	lon_150.append(xy[0])
	lat_150.append(xy[1])

close(f)
dd_box=[array(lon_150).min(), array(lon_150).max(), array(lat_150).min(), array(lat_150).max()]
fat=0.05
box=[dd_box[0]-fat, dd_box[1]+fat, dd_box[2]-fat, dd_box[3]+fat]
f=figure()
mapobj= Basemap(projection='merc',lat_0=(box[2]+box[3])/2.0, lon_0=(box[0]+box[1])/2.0,llcrnrlat=box[2], llcrnrlon=box[0], urcrnrlat=box[3] , urcrnrlon=box[1], resolution='l',area_thresh=1., lat_ts=(box[2]+box[3])/2.0)
longr, latgr=meshgrid(lons,lats)
xx, yy = mapobj(longr, latgr)
mapobj.drawmapboundary()
mapobj.readshapefile('/flurry/home/scollis/shapes/cstntcd_r','coast',drawbounds=True,linewidth=0.5,color='k',antialiased=1,ax=None)
mapobj.contour(xx,yy,angs, levels=[150,30], colors=['r'])
mapobj.drawmeridians(array([130.2, 130.4, 130.6, 130.8,131.0,131.2, 131.4]), labels=[1,0,0,1])
mapobj.drawparallels(array([--12.8, -12.6, -12.4, -12.2, -12.0, -11.8, -11.6, -11.4]), labels=[1,0,0,1])
dd_box_xx, dd_box_yy=mapobj(array([dd_box[0], dd_box[1]]),array([dd_box[2], dd_box[3]]))
xy=dd_box_xx[0], dd_box_yy[0]
width, height= dd_box_xx[1]-dd_box_xx[0], dd_box_yy[1]-dd_box_yy[0]
my_patch=matplotlib.patches.Rectangle(xy, width, height, edgecolor='blue', facecolor='white')
ax=gca()
ax.add_patch(my_patch)
radar_xx, radar_yy=mapobj(array([ber_loc[1], gp_loc[1]]), array([ber_loc[0], gp_loc[0]]))
mapobj.plot(radar_xx, radar_yy, 'bo')
ax.text(radar_xx[0]+1000.0, radar_yy[0]-3000.0, 'Berrimah')
ax.text(radar_xx[1]+1000.0, radar_yy[1]-3000.0, 'Gunn Point')
savefig('/flurry/home/scollis/results/area_vis.png')
close(f)

# 需要导入模块: from matplotlib.toolkits.basemap import Basemap [as 别名]
# 或者: from matplotlib.toolkits.basemap.Basemap import readshapefile [as 别名]
# Tissot's Indicatrix (http://en.wikipedia.org/wiki/Tissot's_Indicatrix). 
# These diagrams illustrate the distortion inherent in all map projections.
# In conformal projections, where angles are conserved around every location, 
# the Tissot's indicatrix are all circles, with varying sizes. In equal-area 
# projections, where area proportions between objects are conserved, the 
# Tissot's indicatrix have all unit area, although their shapes and 
# orientations vary with location.

# adapted from http://www.perrygeo.net/wordpress/?p=4

# create new figure
fig=p.figure()
m = Basemap(llcrnrlon=-180,llcrnrlat=-80,urcrnrlon=180,urcrnrlat=80,
            projection='cyl')
shp_info = m.readshapefile('tissot','tissot',drawbounds=True)
ax = p.gca()
for nshape,seg in enumerate(m.tissot):
    poly = Polygon(seg,facecolor='green',zorder=10)
    ax.add_patch(poly)
# draw meridians and parallels.
m.drawparallels(p.arange(-90,91,30),labels=[1,0,0,0])
m.drawmeridians(p.arange(-180,180,60),labels=[0,0,0,1])
m.drawcoastlines()
m.fillcontinents()
p.title('Tissot Diagram - Cylindrical Equal Area')
print 'plot Cylindrical Equidistant Equal Area Tissot diagram ...'

# create new figure
fig=p.figure()
m = Basemap(llcrnrlon=-180,llcrnrlat=-70,urcrnrlon=180,urcrnrlat=70,