Python Basemap.drawgreatcircle方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def map_trace(tr,ax='None',showpath=True,showplot=True,Lat_0=0.0,Lon_0=60.0):
   from mpl_toolkits.basemap import Basemap

   if ax == 'None':
       m = Basemap(projection='ortho',lat_0=Lat_0,lon_0=Lon_0,resolution='l')
       m.drawmapboundary()
       m.drawcoastlines()
       m.fillcontinents(color='gray',lake_color='white')
   else:
      m = ax

   x1,y1 = m(tr.stats.sac['evlo'],tr.stats.sac['evla'])
   x2,y2 = m(tr.stats.sac['stlo'],tr.stats.sac['stla'])
   m.scatter(x1,y1,s=200.0,marker='*',facecolors='y',edgecolors='k',zorder=99)
   m.scatter(x2,y2,s=20.0,marker='^',color='b',zorder=99)

   if showpath == True:
      m.drawgreatcircle(tr.stats.sac['evlo'],tr.stats.sac['evla'],
                        tr.stats.sac['stlo'],tr.stats.sac['stla'],
                        linewidth=1,color='k',alpha=0.5)

   if showplot == True:
      plt.show()
   else:
      return m

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
 def draw(self, pop):
     fig=plt.figure()
     ax=fig.add_axes([0.1,0.1,0.8,0.8])
     m = Basemap(llcrnrlon=-125.,llcrnrlat=25.,urcrnrlon=-65.,urcrnrlat=52.,
                 rsphere=(6378137.00,6356752.3142),
                 resolution='l',projection='merc',
                 lat_0=40.,lon_0=-20.,lat_ts=20.)
     l = pop[0]          
     for i in range(len(l.sol)):
         lat1 = l.sol[i].lat
         lon1 = l.sol[i].lon
         m.drawgreatcircle(lon1,lat1,lon1,lat1, linewidth=4, color = 'r')
         if i == len(l.sol) - 1:
             lat2 = l.sol[0].lat
             lon2 = l.sol[0].lon
         else:
             lat2 = l.sol[i+1].lat
             lon2 = l.sol[i+1].lon
         m.drawgreatcircle(lon1,lat1,lon2,lat2, color = 'b')
     m.drawcoastlines()
     m.drawstates()
     m.drawcountries()
     m.fillcontinents()
     ax.set_title('GREEDY')
     plt.show()            

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
  def plot_eq(self,ax='None',showpath=True,showplot=True,lon_0=0.0,lat_0=0.0):
     from mpl_toolkits.basemap import Basemap
     
     if ax == 'None':
         #m = Basemap(projection='hammer',lon_0=self.ry)
         m = Basemap(projection='ortho',lat_0=lat_0,lon_0=lon_0,resolution='l')
         m.drawmapboundary()
         m.drawcoastlines()
         m.fillcontinents(color='gray',lake_color='white')
     else:
         m = ax

     x1,y1 = m(self.sy,90.0-self.sx)
     x2,y2 = m(self.ry,90.0-self.rx)
     m.scatter(x1,y1,s=200.0,marker='*',facecolors='y',edgecolors='k',zorder=99)
     m.scatter(x2,y2,s=20.0,marker='^',color='b',zorder=99)
     
     if showpath == True:
        lon_s = self.sy
        lat_s = 90.0-self.sx
        lon_r = self.ry
        lat_r = 90.0-self.rx
        print "lon_s,lat_s,lon_r,lat_r", lon_s, lat_s, lon_r, lat_r
        m.drawgreatcircle(lon_s,lat_s,lon_r,lat_r,linewidth=1,color='k',alpha=0.5)

     if showplot == True:
        plt.show()
     else:
        return m

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def plot_lat_long(coords):
  
  starting_lat = coords[0][1].lat
  starting_lon = coords[0][1].long
  print "starting lat: %s\nstarting lon: %s\n" % (starting_lat, starting_lon)
  for coord in coords[1:]:
    lat = coord[1].lat
    lon = coord[1].long
    print "lat: %s\nlon: %s\n" % (lat, lon)
    m = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80,\
              llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')
    m.drawcoastlines()
    m.fillcontinents(color='coral',lake_color='aqua')

# draw parallels and meridians.
    m.drawparallels(np.arange(-90.,91.,30.))
    m.drawmeridians(np.arange(-180.,181.,60.))
    xpt, ypt =  m(lon, lat)
    m.plot(xpt, ypt, 'bo')
    lonlat = 51.53; lonlon = 0.08
    m.drawgreatcircle(starting_lon,starting_lat,lon,lat,linewidth=2,color='b')
    m.drawmapboundary(fill_color='aqua')
    plt.title("Blorp")
    plt.show()
    starting_lat = lat
    starting_lon = lon

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def general(input_file="input.csv"):
    # create arrays to hold
    lats = []
    longs = []
    weights = []
    names = []

    with open(input_file) as input_csv:
        csv_reader = csv.reader(input_csv)
        #next(csv_reader, None)
        for row in csv_reader:
            #print (float(row[0]))
            
            x, y, w, n = float(row[0]),  float(row[1]), row[2], row[3]
            lats.append(x)
            longs.append(y)
            weights.append(w)
            names.append(n)


    # create background map
    minLat=min(lats)
    maxLat=max(lats)
    minLong = min(longs)
    maxLong = max(longs)    
    if minLat -10 >-90:
        minLat = minLat-10
    if maxLat +10 <90:
        maxLat = maxLat +10
    if minLong -10 >-180:
        minLong = minLong-10
    if maxLong +10<180:
        maxLong = maxLong+10
    m= Basemap(llcrnrlon=minLong, llcrnrlat=minLat, urcrnrlon=maxLong,
                               urcrnrlat=maxLat, projection='merc', lat_1=33, lat_2=45,
                               lon_0=-95, resolution='c', area_thresh=10000)
    m.drawcoastlines()
    m.drawmapboundary(fill_color='aqua')
    m.drawcountries()
    m.fillcontinents(color='green', lake_color='aqua')
    m.drawstates()
    ypt, xpt = m(longs, lats)
    m.plot(ypt, xpt, 'bo', markersize=10)
    for y, x, name in zip(ypt,xpt,names):
        plt.text(y-1000,x+1000,name)

    for i in range(0,len(ypt)-1):
        m.drawgreatcircle(longs[i], lats[i], longs[i+1], lats[i+1], linewidth=1, color='r')
    
    for i in range(0,len(weights)-1):
        plt.text(((ypt[i]+ypt[i+1])/2),((xpt[i]+xpt[i+1])/2),weights[i+1])
    
    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
class GCMapper:
    def __init__(self):
        # create new figure, axes instances.
        self.fig = plt.figure()
        self.ax  = self.fig.add_axes([0.1,0.1,0.8,0.8], axisbg = '#a5bfdd')
        
        # setup mercator map projection.
        self.map = Basemap(llcrnrlon = -130., llcrnrlat = 1.,
                    urcrnrlon = 40.,   urcrnrlat = 70.,
                    rsphere = (6378137.00,6356752.3142),
                    resolution = 'c',
                    projection = 'merc',
                    lat_0 = 40.,lon_0 = -20.,lat_ts = 20.)
        
    def plot(self, element):
        if isinstance(element, Segment):
            self.plot_segment(element)
        if isinstance(element, Route):
            self.plot_route(element)
        if isinstance(element, Point):
            self.plot_point(element)
        raise Exception('Incompatible type %s of element to plot' % type(element))
        
    def plot_route(self, route, color = 'b', linestyle = '-'):
        #just to be suer
        route.init_segments()
        for segment in route.segments:
            self.plot_segment(segment, color = color, linestyle = linestyle)
    
    def plot_point(self, point, formatting = 'go'):
        x, y = self.map(point.lon, point.lat)
        self.map.plot(x, y, formatting, ms = 8)
        
    def plot_segment(self, segment, color = 'b', linestyle = '-'):
        point1 = segment.start
        point2 = segment.end
        self.map.drawgreatcircle(
            point1.lon, point1.lat,
            point2.lon, point2.lat,
            linestyle = linestyle,
            linewidth = 1.5,
            color = color)

    def render(self):
        self.map.drawcoastlines(color='#8f8457')
        self.map.fillcontinents(color='#f5f0db')
        self.map.drawcountries(color='#a9a06d')
        self.map.drawparallels(np.arange(10,90,20), labels = [1,1,0,1])
        self.map.drawmeridians(np.arange(-180,180,30), labels = [1,1,0,1])
        self.ax.set_title('Flights')
        plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
   def find_pierce_coor(self,plot='False'):
      import geopy
      from geopy.distance import VincentyDistance

      '''
      given an instance of the receiver function class this function
      returns latitude and longitude of all receiver side pierce points
      of Pds in a given depth range (the default range is 50 - 800 km)
      NOTE:
      be careful that ses3d typically uses colatitude, while
      this function returns latitude '''

      depth_range = np.arange(50,800,5)        #set range of pierce points

      #geodetic info
      bearing     = self.az
      lon_s = self.ses3d_seismogram.sy
      lat_s = 90.0-self.ses3d_seismogram.sx
      lon_r = self.ses3d_seismogram.ry
      lat_r = 90.0-self.ses3d_seismogram.rx
      origin      = geopy.Point(lat_s, lon_s)

      #find how far away the pierce point is
      model  = TauPyModel(model='pyrolite_5km')

      for i in range(0,len(depth_range)):
         phase = 'P'+str(depth_range[i])+'s'
         pierce = model.get_pierce_points(self.eq_depth,self.delta_deg,phase_list=[phase])
         points = pierce[0].pierce
         for j in range(0,len(points)):
            if points[j]['depth'] == depth_range[i] and points[j]['dist']*(180.0/np.pi) > 25.0:
               prc_dist = points[j]['dist']*(180.0/np.pi)
               d_km = prc_dist * ((2*np.pi*6371.0/360.0))
               destination = VincentyDistance(kilometers=d_km).destination(origin,bearing)
               lat = destination[0]
               lon = destination[1]
               value = 0
               row = {'depth':depth_range[i],'dist':prc_dist,'lat':lat,'lon':lon,'value':value}
               self.pierce_dict.append(row)

      if plot=='True':
         m = Basemap(projection='hammer',lon_0=0)
         m.drawmapboundary()
         m.drawcoastlines()
         m.drawgreatcircle(lon_s,lat_s,lon_r,lat_r,linewidth=1,color='b',alpha=0.5)

         for i in range(len(self.pierce_dict)):
            x,y = m(self.pierce_dict[i]['lon'],self.pierce_dict[i]['lat'])
            m.scatter(x,y,5,marker='o',color='r')
         plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
    def _dump_graph(self, title=None, lon_0=0, lat_0=90,
                    projection='vandg', func=lambda x: len(x._edges)):
        from mpl_toolkits.basemap import Basemap
        from matplotlib import pyplot as plt
        fig = plt.figure()
        m = Basemap()

        counts = {}
        for node in self._nodes:
            count = func(node)
            counts.setdefault(count, [])
            counts[count].append(list(node._point))

        minx = np.inf
        miny = np.inf
        maxx = -np.inf
        maxy = -np.inf
        for k, v in counts.items():
            v = np.array(v)
            ra, dec = vector.vector_to_radec(v[:, 0], v[:, 1], v[:, 2])
            x, y = m(ra, dec)
            m.plot(x, y, 'o', label=str(k))
            for x0 in x:
                minx = min(x0, minx)
                maxx = max(x0, maxx)
            for y0 in y:
                miny = min(y0, miny)
                maxy = max(y0, maxy)

        for edge in list(self._edges):
            A, B = [x._point for x in edge._nodes]
            r0, d0 = vector.vector_to_radec(A[0], A[1], A[2])
            r1, d1 = vector.vector_to_radec(B[0], B[1], B[2])
            m.drawgreatcircle(r0, d0, r1, d1, color='blue')

        plt.xlim(minx, maxx)
        plt.ylim(miny, maxy)
        if title:
            plt.title("%s, %d v, %d e" % (
                title, len(self._nodes), len(self._edges)))
        plt.legend()
        plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
    def GDELT_interactions_maplot(self,counts_int):

        max_val = np.log10(counts_int.max())

        def get_alpha(interaction_counts):
            '''
            Convert a count to an alpha val.
            Log-scaled
            '''
            scale = np.log10(interaction_counts)
            return (scale/max_val) * 0.25

        # Draw the basemap like before
        plt.figure(figsize=(12,12))
        event_map = Basemap(projection='merc',
                            resolution='l', area_thresh=1000.0,  # Low resolution
                            lat_0=15, lon_0=30,  # Map center
                            llcrnrlon=10, llcrnrlat=1,  # Lower left corner
                            urcrnrlon=50, urcrnrlat=30)  # Upper right corner
        # Draw important features
        event_map.drawcoastlines()
        event_map.drawcountries()
        event_map.fillcontinents(color='0.8')
        event_map.drawmapboundary()

        # Draw the lines on the map:
        for arc, count in interaction_counts.iteritems():
            point1, point2 = arc
            y1, x1 = point1
            y2, x2 = point2

            # Only plot lines where both points are on our map:
            if ((x1 > 10 and x1 < 100 and y1 > 20 and y1 < 70) and
                (x2 > 10 and x2 < 100 and y2 > 20 and y2 < 70)):

                line_alpha = get_alpha(count)

                # Draw the great circle line
                event_map.drawgreatcircle(x1, y1, x2, y2, linewidth=2,color='r', alpha=line_alpha)
        plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def translation_pathway(N, vt=2., Ric=1221):

    fig, ax = plt.subplots(1)
    fig2, ax2 = plt.subplots(1)
    # use low resolution coastlines.
    map = Basemap(projection='moll',lat_0=0,lon_0=0,resolution='l')
    # draw coastlines, country boundaries, fill continents.
    map.drawcoastlines(linewidth=0.25)
    map.drawcountries(linewidth=0.25)
    map.fillcontinents(color='coral',lake_color='aqua')
    # draw the edge of the map projection region (the projection limb)
    map.drawmapboundary(fill_color='aqua')
    # draw lat/lon grid lines every 30 degrees.
    map.drawmeridians(np.arange(0,360,30))
    map.drawparallels(np.arange(-90,90,30))

    for i in range(N):
        
        Point = seismo.random_points(seismo="surface")
        pos_1, pos_2 = seismo.bottom_point_zeta(Point, Ric=1221)

        x_, y_, z_ = seismo.from_seismo_to_cartesian(Ric-Point['depth'],
                                                  Point['longitude'], Point['latitude'])
        Age = geodyn.translation(x_/Ric, y_/Ric, z_/Ric, vt, 1.)
        age_x, age_y = map(Point['longitude'], Point['latitude'])
        map.scatter(age_x, age_y, s=Age*100, zorder=10)

        #print Age, Ric-Point['depth'], Point['longitude'], Point['longitude']
        Nd = 20
        x, y, z, d, dx = seismo.raypath_straight(pos_1, pos_2, Nd, coordinate_type="spherical")
        r, theta, phi = seismo.from_cartesian_to_seismo(x, y, z)
        
        A = geodyn.translation(x/Ric, y/Ric, z/Ric, vt, 1.)
        Age_average = np.sum(A)/Nd
        
        map.scatter(age_x, age_y, s=Age*100, zorder=20)
        ax.scatter(Point['longitude'], Age_average)
        map.drawgreatcircle(pos_1['longitude'], pos_1['latitude'], pos_2['longitude'], pos_2['latitude'])
        for i, xtheta in enumerate(theta[0:-1]):
            map.drawgreatcircle(phi[i], theta[i], phi[i+1], theta[i+1])

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def example_2():
    
    # retrieve data
    dir  = "/homespace/gaubert/viirs/real-data"
    #geo_file     = h5py.File("%s/%s" %(dir,"GMODO_npp_d20120224_t1821456_e1823098_b01694_c20120306214722023793_cspp_dev.h5"))
    #geo_file     = h5py.File("%s/%s" %(dir,"GMODO_npp_d20120224_t1823110_e1824352_b01694_c20120306215057805679_cspp_dev.h5"))
    #geo_file     = h5py.File("%s/%s" %(dir,"GMODO_npp_d20120224_t1826018_e1827260_b01694_c20120306215454419051_cspp_dev.h5"))
    #geo_file     = h5py.File("%s/%s" %(dir,"GMODO_npp_d20120224_t1827272_e1828514_b01694_c20120306215852012949_cspp_dev.h5"))
    geo_file      = h5py.File("/homespace/gaubert/GMTCO_npp_d20120224_t1100479_e1102121_b01689_c20120224172231282331_noaa_ops.h5")
    
    
    #lats = geo_file['All_Data']['VIIRS-MOD-GEO_All']['Latitude'][:]
    #lons = geo_file['All_Data']['VIIRS-MOD-GEO_All']['Longitude'][:]
    lats = geo_file['All_Data']['VIIRS-MOD-GEO-TC_All']['Latitude'][:]
    lons = geo_file['All_Data']['VIIRS-MOD-GEO-TC_All']['Longitude'][:]
    
    line_len = len(lats[0])
    col_len  = len(lats)
    
    print("line len %d, col len %d" % (line_len, col_len))
    print("Upper left corner point: (%f,%f)\n" % (lats[0][0], lons[0][0] ))
    print("Lower right corner point: (%f,%f)\n" % (lats[col_len-1][line_len-1], lons[col_len-1][line_len-1]))
    
    
    # draw map with markers for float locations
    #m = Basemap(projection='hammer',lon_0=180)
    lon_ref = lons[(col_len-1)/2][(line_len-1)/2]
    #lat_ref = 10
    lat_ref = lats[(col_len-1)/2][(line_len-1)/2]
    #m = Basemap(projection='ortho',lat_0=lat_ref,lon_0=lon_ref,resolution='l')
    m = Basemap(projection='nsper',lat_0=lat_ref,lon_0=lon_ref,satellite_height=2000*1000,resolution='l')
    
    #x, y = m(lons[0:10],lats[0:10])
    x,y  = m(lons,lats)
    
    
    m.drawcoastlines()
    
    m.drawmapboundary(fill_color='#99ffff')
    #m.fillcontinents(color='#cc9966',lake_color='#99ffff')
    #m.scatter(x,y,s = 1 ,color='k')
    
    m.drawgreatcircle(lons[0][0],lats[0][0],lons[0][-1],lats[0][-1],linewidth=1,color='b')
    m.drawgreatcircle(lons[0][0],lats[0][0],lons[col_len-1][0],lats[col_len-1][0],linewidth=1,color='b')
    m.drawgreatcircle(lons[col_len-1][0],lats[col_len-1][0],lons[col_len-1][line_len-1],lats[col_len-1][line_len-1],linewidth=1,color='b')
    m.drawgreatcircle(lons[0][line_len-1],lats[0][line_len-1],lons[col_len-1][line_len-1],lats[col_len-1][line_len-1],linewidth=1,color='b')
    
    plt.title('Location of VIIRS granule ')
    plt.savefig('/tmp/plot-gran3.png')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def draw(x1,x2,x3,x4):

	map = Basemap(projection='robin', lat_0=0, lon_0=0, resolution='l', area_thresh=1000.0)
	map.drawcoastlines()
	map.drawcountries()
	map.fillcontinents(color='coral', lake_color='aqua')
	map.drawmapboundary(fill_color='aqua')

	map.drawmeridians(np.arange(0,360,30))
	map.drawparallels(np.arange(-90,90,30))

	#initiate ocean and icecap colours
	#map.bluemarble()

	#draw a blob on the map
	x,y = map(x2,x1)
	map.plot(x,y, 'bo', markersize=10)

	x,y = map(x4,x3)
	map.plot(x,y, 'bd', markersize=10)
	
	map.drawgreatcircle(x2,x1,x4,x3,linewidth=2,color='b')

	plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def map_stream(st,showpath=True,showplot=True,Lat_0=0.0,Lon_0=60.0):
   from mpl_toolkits.basemap import Basemap

   fig = plt.figure()
   plt.ion()
   m = Basemap(projection='ortho',lat_0=Lat_0,lon_0=Lon_0,resolution='l')
   m.drawmapboundary()
   m.drawcoastlines()
   m.fillcontinents(color='gray',lake_color='white')
   gc_lines = []

   for tr in st:
      x1,y1 = m(tr.stats.sac['evlo'],tr.stats.sac['evla'])
      x2,y2 = m(tr.stats.sac['stlo'],tr.stats.sac['stla'])
      m.scatter(x1,y1,s=200.0,marker='*',facecolors='y',edgecolors='k',zorder=99)
      station_pt = m.scatter(x2,y2,s=20.0,marker='^',color='b',zorder=99,picker=1)
      station_pt.set_label(tr.stats.station)

      if showpath == True:
         gc_line = m.drawgreatcircle(tr.stats.sac['evlo'],tr.stats.sac['evla'],
                                     tr.stats.sac['stlo'],tr.stats.sac['stla'],
                                     linewidth=1,color='k',alpha=0.5)
         gc_line[0].set_label(tr.stats.station)
         gc_lines.append(gc_line)

   def onpick(event):
      art = event.artist
      picked =  art.get_label()
      print "removing station(s) ", picked
      st_to_remove = st.select(station=picked)
      for killed_tr in st_to_remove:
         st.remove(killed_tr)
      art.remove()

      for gc in gc_lines:
         gc_line_label = gc[0].get_label()
         if gc_line_label == picked:
            gc[0].remove()

      fig.canvas.draw()
   
   fig.canvas.mpl_connect('pick_event', onpick)

   if showplot == True:
      plt.show()
   else:
      return m

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
def save_map(routes, file_name):
    """
    Render flight routes to an image file.

    :param list routes: flight path lines
    :param str file_name: image output file name
    """
    fig = plt.figure(figsize=(7.195, 3.841), dpi=100)
    m = Basemap(projection='cyl', lon_0=0, resolution='c')
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)

    for (colour, alpha, linewidth, lat1, long1, lat2, long2) in sorted(routes):
        """
        Cannot handle situations in which the great circle intersects the 
        edge of the map projection domain, and then re-enters the domain.

        Fix from: http://stackoverflow.com/questions/13888566/
        """
        line, = m.drawgreatcircle(long1, lat1, long2, lat2,
                                  linewidth=linewidth,
                                  color=colour,
                                  alpha=alpha,
                                  solid_capstyle='round')

        p = line.get_path()

        # Find the index which crosses the dateline (the delta is large)
        cut_point = np.where(np.abs(np.diff(p.vertices[:, 0])) > 200)[0]

        if cut_point:
            cut_point = cut_point[0]

            # Create new vertices with a nan in between and set 
            # those as the path's vertices
            new_verts = np.concatenate([p.vertices[:cut_point, :], 
                                        [[np.nan, np.nan]], 
                                        p.vertices[cut_point+1:, :]])
            p.codes = None
            p.vertices = new_verts

    m.warpimage(image="earth_lights_lrg.jpg")
    plt.savefig(file_name, dpi=1000)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import drawgreatcircle [as 别名]
            rsphere=(6378137.00,6356752.3142),\
            resolution='c',area_thresh=10000.,projection='merc',\
            lat_0=40.,lon_0=-20.,lat_ts=20.)
# nylat, nylon are lat/lon of New York
nylat = 40.78
nylon = -73.98
# lonlat, lonlon are lat/lon of London.
lonlat = 51.53
lonlon = 0.08

# find 1000 points along the great circle.
#x,y = m.gcpoints(nylon,nylat,lonlon,lonlat,1000)
# draw the great circle.
#m.plot(x,y,linewidth=2)
# drawgreatcircle performs the previous 2 steps in one call.
m.drawgreatcircle(nylon,nylat,lonlon,lonlat,linewidth=2,color='b')

m.drawcoastlines()
m.fillcontinents()
# draw parallels
circles = np.arange(10,90,20)
m.drawparallels(circles,labels=[1,1,0,1])
# draw meridians
meridians = np.arange(-180,180,30)
m.drawmeridians(meridians,labels=[1,1,0,1])
plt.title('Great Circle from New York to London (Mercator)')
sys.stdout.write('plotting Great Circle from New York to London (Mercator)\n')

# create new figure
fig=plt.figure()
# setup a gnomonic projection.