Python mpl_toolkits.basemap 模块，Basemap() 实例源码

def plot_grid2D(lons, lats, tec_grid2D, datetime, title_label = ''):

    LATS, LONS = np.meshgrid(lats, lons)

    m = Basemap(llcrnrlon=-180,
                llcrnrlat=-55,
                urcrnrlon=180,
                urcrnrlat=75,
                projection='merc',
                area_thresh=1000,
                resolution='i')

    m.drawstates()
    m.drawcountries()
    m.drawcoastlines()

    parallels = np.arange(-90,90,20)
    m.drawparallels(parallels,labels=[True,False,False,True])
    meridians = np.arange(0,360,40)
    m.drawmeridians(meridians,labels=[True,False,False,True])

    m.scatter(LONS, LATS, c=tec_grid2D, latlon = True, linewidths=0, s=5)
    m.colorbar()

    plt.title('%s\n%s' % (title_label, datetime.isoformat(' ')))

def prepare_basemap(min_lat, min_lon, max_lat, max_lon, delta_lat, delta_lon):
    """

    :param min_lat: float 
    :param min_lon: float
    :param max_lat: float
    :param max_lon: float
    :param delta_lat: float
    :param delta_lon: float
    :return: A Basemap instance
    """
    m = Basemap(projection='cyl', llcrnrlat=min_lat, urcrnrlat=max_lat, urcrnrlon=max_lon,
                llcrnrlon=min_lon, resolution='l')
    m.drawcountries()
    m.drawcoastlines()
    m.drawparallels(np.arange(np.around(min_lat, 1), np.around(max_lat, 1), -round(delta_lat) / 5.0),
                    labels=[1, 0, 0, 0])
    m.drawmeridians(np.arange(np.around(min_lon, 1), np.around(max_lon, 1), -round(delta_lon) / 5.0),
                    labels=[0, 0, 0, 1])
    return m

def make_stations(network,latmin,latmax,lonmin,lonmax,dlat,dlon,plot=True,**kwargs):
    head = kwargs.get('head','RM')
    elev = kwargs.get('elev',0)
    dep  = kwargs.get('dep',0)

    lats = np.arange(latmin,latmax+dlat,dlat)
    lons = np.arange(lonmin,lonmax+dlon,dlon)
    lats_pts,lons_pts = np.meshgrid(lats,lons)

    if plot:
        m = Basemap(projection='hammer',lon_0=0,resolution='l')
        m.drawcoastlines()
        m.drawmeridians(np.arange(0,351,10))
        m.drawparallels(np.arange(-80,81,10))
        lons_pts2,lats_pts2 = m(lons_pts,lats_pts)
        m.scatter(lons_pts2,lats_pts2)
        plt.show()

    f = open('STATIONS','w')
    print 'total number of stations : ',len(lats_pts)
    st = 1
    for i in range(0,len(lats)):
        for j in range(0,len(lons)):
            f.write('{}{:04d} {} {:5.2f} {:5.2f} {:5.2f} {:5.2f}'.format(head,st,network,lats[i],lons[j],elev,dep)+'\n')
            st += 1

def plot_geo_config(stations,events):
   '''
   plot geometrical configuration of tomography experiment
   '''
   m = Basemap(projection='hammer',lon_0=0,resolution='l')
   m.drawcoastlines()
   m.drawmeridians(np.arange(0,361,30))
   m.drawparallels(np.arange(-90,91,30))
   station_lons = stations[0,:]
   station_lats = stations[1,:]
   x,y = m(station_lons,station_lats)
   m.scatter(x,y,s=50,marker='^',c='r')

   for event in events:
      lon = event[1]
      lat = event[2]
      x,y = m(lon,lat)
      m.scatter(x,y,s=100,marker='*',c='y')

   plt.savefig('geo_config.pdf',format='pdf')

def plot_pierce_points(self,depth=410.0,ax='None'):
      '''
      Plots pierce points for a given depth.  If an axes object is supplied
      as an argument it will use it for the plot.  Otherwise, a new axes 
      object is created.
      '''

      if ax == 'None':
         m = Basemap(llcrnrlon=0.0,llcrnrlat=-35.0,urcrnrlon=120.0,urcrnrlat=35.0)
      else:
         m = ax

      m.drawmapboundary()
      m.drawcoastlines()
      found_points = False

      for i in range(0,len(self.pierce_dict)):
         if self.pierce_dict[i]['depth'] == depth:
            x,y = m(self.pierce_dict[i]['lon'],self.pierce_dict[i]['lat'])
            found_points = True

      if found_points == True:
         m.scatter(x,y,50,marker='+',color='r')
      else:
         print "no pierce points found for the given depth"

def plot_contiguous_US_tweets(lon, lat, file_path):
    '''
    INPUT: List of longitudes (lon), list of latitudes (lat), file path to save the plot (file_path) 
    OUTPUT: Plot of tweets in the contiguous US.
    '''

    map = Basemap(projection='merc',
                  resolution = 'h', 
                  area_thresh = 10000,
                  llcrnrlon=-140.25,    # lower left corner longitude of contiguous US
                  llcrnrlat=5.0,        # lower left corner latitude of contiguous US
                  urcrnrlon=-56.25,     # upper right corner longitude of contiguous US
                  urcrnrlat=54.75)      # upper right corner latitude of contiguous US

    x,y = map(lon, lat)

    map.plot(x, y, 'bo', markersize=2, alpha=.3)
    map.drawcoastlines()
    map.drawstates()
    map.drawcountries()
    map.fillcontinents(color = '#DAF7A6', lake_color='#a7cdf2')
    map.drawmapboundary(fill_color='#a7cdf2')
    plt.gcf().set_size_inches(15,15)
    plt.savefig(file_path, format='png', dpi=1000)

def main2():
    from mpl_toolkits.basemap import Basemap
    import numpy as np
    import matplotlib.pyplot as plt
    # llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
    # are the lat/lon values of the lower left and upper right corners
    # of the map.
    # resolution = 'i' means use intermediate resolution coastlines.
    # lon_0, lat_0 are the central longitude and latitude of the projection.
    m = Basemap(llcrnrlon=9.5, llcrnrlat=63.2,
                urcrnrlon=10.9, urcrnrlat=64.,
                resolution='i', projection='tmerc',
                lon_0=10.7, lat_0=63.4)
    # can get the identical map this way (by specifying width and
    # height instead of lat/lon corners)
    # m = Basemap(width=894887,height=1116766,\
    #            resolution='i',projection='tmerc',lon_0=-4.36,lat_0=54.7)
    m.drawcoastlines()
    m.fillcontinents(color='coral', lake_color='aqua')
    # m.drawparallels(np.arange(-40, 61., 2.))
    # m.drawmeridians(np.arange(-20., 21., 2.))
    m.drawmapboundary(fill_color='aqua')
    plt.title("Transverse Mercator Projection")
    plt.show()

def _initialize_background_map(self):  # calculate the correct aspect ratio by drawing a map of the substrate
        fig = plt.figure()
        ax = plt.subplot(111)
        nodes_to_map_coordinates = [MapCoordinate(lat=data['Latitude'], lon=data['Longitude'])
                                    for n, data in self._substrate.node.items()]
        lower_left, upper_right = ExtendedGraphVisualizer._get_corner_map_coordinates(nodes_to_map_coordinates)
        self._background_map = Basemap(resolution="c",
                                       projection='merc',
                                       anchor="W",
                                       ax=ax,
                                       llcrnrlat=lower_left.lat,
                                       urcrnrlat=upper_right.lat,
                                       llcrnrlon=lower_left.lon,
                                       urcrnrlon=upper_right.lon,
                                       fix_aspect=True)
        self._background_map.fillcontinents(ax=ax)
        (xmin, xmax), (ymin, ymax) = ax.get_xlim(), ax.get_ylim()
        dx, dy = float(abs(xmax - xmin)), float(abs(ymax - ymin))
        self._map_width = (dx / dy) * MAP_HEIGHT
        plt.close(fig)

def plot_map():
    lon=grd.lon; lat=grd.lat
    lat_ts=((lat[0,0]+lat[-1,0])*0.5)
    cbox=[lon[0, 0],lat[0, 0],lon[-1, -1],lat[-1, -1]]
    map = Basemap(projection='merc', resolution='l',
                    lat_ts=lat_ts,
                    llcrnrlon=cbox[0], llcrnrlat=cbox[1],
                    urcrnrlon=cbox[2], urcrnrlat=cbox[3])
    map.drawmapboundary(fill_color='white')
    map.drawcoastlines()
    map.fillcontinents(color='#cc9966', lake_color='#99ffff')
    lon_tcks = np.arange(-100.,30.,20.)
    lat_tcks = np.arange(0.,70.,10.)
    map.drawmeridians(lon_tcks, labels=[True, False, False, True])
    map.drawparallels(lat_tcks, labels=[False, True, True, False])
    # get projected coordinates
    x,y=map(lon,lat)
    return map, x, y

def basemap_raster_mercator(lon, lat, grid, cmin, cmax, cmap_name):

    # longitude/latitude extent
    lons = (np.amin(lon), np.amax(lon))
    lats = (np.amin(lat), np.amax(lat))

    # construct spherical mercator projection for region of interest
    m = Basemap(projection='merc',llcrnrlat=lats[0], urcrnrlat=lats[1],
                llcrnrlon=lons[0],urcrnrlon=lons[1])

    #vmin,vmax = np.nanmin(grid),np.nanmax(grid)
    masked_grid = np.ma.array(grid,mask=np.isnan(grid))
    fig = plt.figure(frameon=False,figsize=(12,8),dpi=72)
    plt.axis('off')
    cmap = mpl.cm.get_cmap(cmap_name)
    m.pcolormesh(lon,lat,masked_grid,latlon=True,cmap=cmap,vmin=cmin,vmax=cmax)

    str_io = StringIO.StringIO()
    plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True)
    plt.close()

    numpy_bounds = [ (lons[0],lats[0]),(lons[1],lats[0]),(lons[1],lats[1]),(lons[0],lats[1]) ]
    float_bounds = [ (float(x), float(y)) for x,y in numpy_bounds ]
    return str_io.getvalue(), float_bounds

def basemap_barbs_mercator(u,v,lat,lon):

    # lon/lat extents
    lons = (np.amin(lon), np.amax(lon))
    lats = (np.amin(lat), np.amax(lat))

    # construct spherical mercator projection for region of interest
    m = Basemap(projection='merc',llcrnrlat=lats[0], urcrnrlat=lats[1],
                llcrnrlon=lons[0],urcrnrlon=lons[1])

    #vmin,vmax = np.nanmin(grid),np.nanmax(grid)
    fig = plt.figure(frameon=False,figsize=(12,8),dpi=72*4)
    plt.axis('off')
    m.quiver(lon,lat,u,v,latlon=True)

    str_io = StringIO.StringIO()
    plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True)
    plt.close()

    numpy_bounds = [ (lons[0],lats[0]),(lons[1],lats[0]),(lons[1],lats[1]),(lons[0],lats[1]) ]
    float_bounds = [ (float(x), float(y)) for x,y in numpy_bounds ]
    return str_io.getvalue(), float_bounds

def __init__(self,ASN):
        self.filename = str(ASN) + 'network graph'
        self.alias_filename = str(ASN) + 'alias_candidates.txt'
        self.ISP_Network = networkx.Graph()
        self.files_read = 0
        self.border_points = set()
        if os.path.isfile(self.filename):
            load_file = shelve.open(self.filename, 'r')
            self.files_read = load_file['files_read']
            self.ISP_Network = load_file['ISP_Network']
            self.border_points = load_file['border_points']
        self.ASN = ASN
        plt.ion()
        country = raw_input('enter country to focus on map [Israel/Usa/Australia/other]: ')
        if country == 'Israel' or country == 'israel' or country == 'ISRAEL':
            self.wmap = Basemap(projection='aeqd', lat_0 = 31.4, lon_0 = 35, width = 200000, height = 450000, resolution = 'i')
        elif country == 'USA' or country == 'usa':
            self.wmap = Basemap(projection='aeqd', lat_0 = 40, lon_0 = -98, width = 4500000, height = 2700000, resolution = 'i')
        elif country == 'Australia' or 'australia' or 'AUSTRALIA':
            self.wmap = Basemap(projection='aeqd', lat_0 = -23.07, lon_0 = 132.08, width = 4500000, height = 3500000, resolution = 'i')
        else:
            self.wmap = Basemap(projection='cyl', resolution = 'c')
        plt.hold(False)

def make_a_basemap(ax, ll_cnr, ru_cnr):
        bm = Basemap(
            projection=Choroplether.EPSG_4283_APPROXIMATION['projection'],
            ellps=Choroplether.EPSG_4283_APPROXIMATION['ellps'],
            lon_0=(ru_cnr[0] - ll_cnr[0]) / 2,
            lat_0=(ru_cnr[1] - ll_cnr[1]) / 2,
            llcrnrlon=ll_cnr[0],
            llcrnrlat=ll_cnr[1],
            urcrnrlon=ru_cnr[0],
            urcrnrlat=ru_cnr[1],
            lat_ts=0,
            resolution='i',
            suppress_ticks=True,
            ax=ax
        )
        return bm

def get_basemap(shapefile):
    shp = fiona.open(shapefile+'.shp')
    bds = shp.bounds
    shp.close()
    extra = 0.01
    ll = (bds[0], bds[1])
    ur = (bds[2], bds[3])
    coords = list(chain(ll, ur))
    w, h = coords[2] - coords[0], coords[3] - coords[1]
    m = Basemap(projection='tmerc',
                lon_0= -118.2437,
                lat_0= 34.0522,
                ellps = 'WGS84',
                llcrnrlon=coords[0],
            #     llcrnrlat=coords[1] - extra + 0.01 * h,
                llcrnrlat=33.6,
                urcrnrlon=coords[2],
                urcrnrlat=coords[3],
                lat_ts=0,
                resolution='i',
                suppress_ticks=True)
    m.readshapefile(shapefile, 'LA', color='black', zorder=2)
    return m, coords

# Convenience functions for working with colour ramps and bars

def main():
    bm = Basemap(projection = "rotpole",
                      o_lat_p = 36.0,
                      o_lon_p = 180.0,
                      llcrnrlat = -10.590603,
                      urcrnrlat = 46.591976,
                      llcrnrlon = -139.08585,
                      urcrnrlon = 22.661009,
                      lon_0 = -106.0,
                      rsphere = 6370000,
                      resolution = 'l')

    fig = plt.figure(figsize=(8,8))
    ax = fig.add_axes([0.1,0.1,0.8,0.8])

    bm.drawcoastlines(linewidth=.5)

    print bm.proj4string

    plt.savefig("basemap_map.png")
    plt.close(fig)

def _basemap(self, geobounds, **kwargs):
        if not basemap_enabled():
            return None

        local_kwargs = dict(projection = "lcc",
                            lon_0 = self.stand_lon,
                            lat_0 = self.moad_cen_lat,
                            lat_1 = self.truelat1,
                            lat_2 = self.truelat2,
                            llcrnrlat = geobounds.bottom_left.lat,
                            urcrnrlat = geobounds.top_right.lat,
                            llcrnrlon = geobounds.bottom_left.lon,
                            urcrnrlon = geobounds.top_right.lon,
                            rsphere = Constants.WRF_EARTH_RADIUS,
                            resolution = 'l')
        local_kwargs.update(kwargs)

        _basemap = Basemap(**local_kwargs)

        return _basemap

def _basemap(self, geobounds, **kwargs):
        if not basemap_enabled():
            return None

        local_kwargs = dict(projection = "merc",
                            lon_0 = self.stand_lon,
                            lat_0 = self.moad_cen_lat,
                            lat_ts = self._lat_ts,
                            llcrnrlat = geobounds.bottom_left.lat,
                            urcrnrlat = geobounds.top_right.lat,
                            llcrnrlon = geobounds.bottom_left.lon,
                            urcrnrlon = geobounds.top_right.lon,
                            rsphere = Constants.WRF_EARTH_RADIUS,
                            resolution = "l")
        local_kwargs.update(kwargs)

        _basemap = Basemap(**local_kwargs)

        return _basemap

def _basemap(self, **kwargs):
        if not basemap_enabled():
            return None

        local_kwargs = dict(projection = "stere",
                            lon_0 = self.stand_lon,
                            lat_0 = self._hemi,
                            lat_ts = self._lat_ts,
                            llcrnrlat = geobounds.bottom_left.lat,
                            urcrnrlat = geobounds.top_right.lat,
                            llcrnrlon = geobounds.bottom_left.lon,
                            urcrnrlon = geobounds.top_right.lon,
                            rsphere = Constants.WRF_EARTH_RADIUS,
                            resolution = "l")
        local_kwargs.update(kwargs)

        _basemap = Basemap(**local_kwargs)

        return _basemap

def _basemap(self, geobounds, **kwargs):
        if not basemap_enabled():
            return None

        local_kwargs = dict(projection = "rotpole",
                            o_lat_p = self._bm_cart_pole_lat,
                            o_lon_p = self.pole_lon,
                            llcrnrlat = geobounds.bottom_left.lat,
                            urcrnrlat = geobounds.top_right.lat,
                            llcrnrlon = geobounds.bottom_left.lon,
                            urcrnrlon = geobounds.top_right.lon,
                            lon_0 = self._bm_lon_0,
                            rsphere = Constants.WRF_EARTH_RADIUS,
                            resolution = "l")

        local_kwargs.update(kwargs)
        _basemap = Basemap(**local_kwargs)

        return _basemap

def create_basemap(self, projection='merc'):
        """Creates a very basic basemap: 'bbox': (x1, x2, y1, y2)

        If you want to use the power of basemap you can write your own basemap
        and pass it to you object
        """
        bm = Basemap(
            llcrnrlat=self.bbox[2], urcrnrlat=self.bbox[3],
            llcrnrlon=self.bbox[0], urcrnrlon=self.bbox[1],
            projection=projection)
        bm.drawcoastlines(linewidth=0)
        return bm

def create_map():
    """
    create the base map for the choropleth.
    """
    logging.debug('create_map() -- Please wait while I create the world.')
    degrees_width = 118.0
    degrees_height = 55.0
    center_lat = 44.5
    center_lon = -110.0
    my_map = Basemap(  # ax=ax,
        projection='merc',  # default is cyl
        ellps='WGS84',
        lat_0=center_lat, lon_0=center_lon,
        llcrnrlat=center_lat - degrees_height / 2.0,
        llcrnrlon=center_lon - degrees_width / 2.0,
        urcrnrlat=center_lat + degrees_height / 2.0,
        urcrnrlon=center_lon + degrees_width / 2.0,
        resolution='i',  # 'c', 'l', 'i', 'h', 'f'
    )
    logging.debug('created map')
    logging.debug('loading shapes...')
    for section_name in CONTEST_SECTIONS.keys():
        # logging.debug('trying to load shape for %s', section_name)
        try:
            my_map.readshapefile('shapes/%s' % section_name, section_name, drawbounds=False)
        except IOError, err:
            logging.error('Could not load shape for %s' % section_name)

    logging.debug('loaded section shapes')
    return my_map

def Heatmap(self, grib_object, file_name):
        """Generate Heatmap with Data from grib object
        grib_object: an object containing raw data to be visualized
        file_name:   a string representing the name of generated picture
        """
        data = grib_object.values
        lat,lon = grib_object.latlons()
        unit = grib_object['units']
        data_type = grib_object['name']
        date = self.GetTime(grib_object)

        fig = plt.figure(figsize=(8,8))
        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])    

        m = Basemap(
                resolution='c', # c, l, i, h, f or None
                projection='cyl',
                lat_0=39.72, lon_0=-86.29,
                llcrnrlon=-87.79, llcrnrlat= 38.22,
                urcrnrlon=-84.79, urcrnrlat=41.22)

        parallels = np.arange(38.22, 41.22, 0.5)
        m.drawparallels(parallels,labels=[False,True,True,False])
        meridians = np.arange(-87.79, -84.79, 0.5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        x,y = m(lon, lat)
        cs = m.pcolormesh(x,y,data,
                        shading='flat',
                        cmap=plt.cm.jet)

        cbar = plt.colorbar(cs,location='bottom', fraction=0.046, pad=0.06)

        # Adjust the position of Unit
        cbar_ax = cbar.ax
        cbar_ax.text(0.0, -1.3, unit, horizontalalignment='left')
        m.readshapefile(self.shapeFile,'areas')
        plt.title(data_type + '   ' + date, fontsize = 'x-large')
        plt.savefig(file_name)
        plt.close(fig)

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

#########################################################################
# Plot map of earthquake and station for all traces in stream
#########################################################################

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

def draw_basemap(resolution='l', ax=None,country_linewidth=0.5, coast_linewidth=
                     1.0, zorder=None,  **kwds):
    if ax is None:
        ax = plt.gca()
    m = Basemap(*(ax.viewLim.min + ax.viewLim.max), resolution=resolution, ax=ax, **kwds)
    m.drawcoastlines(linewidth=coast_linewidth, zorder=zorder)
    m.drawcountries(linewidth=country_linewidth, zorder=zorder)
    return m

def _create_map(self, coordinates, labels, n_clusters):
        """
        Graphs GPS coordinates on map (used ti verify clustering algo)
        """
        map = Basemap(llcrnrlat=22, llcrnrlon=-119, urcrnrlat=49, urcrnrlon=-64,
        projection='lcc',lat_1=33,lat_2=45,lon_0=-95)

        base_color = 'white'
        border_color = 'lightgray'
        boundary_color = 'gray'

        map.fillcontinents(color=base_color, lake_color=border_color)
        map.drawstates(color=border_color)
        map.drawcoastlines(color=border_color)
        map.drawcountries(color=border_color)
        map.drawmapboundary(color=boundary_color)


        lat, lon = zip(*coordinates)
        map.scatter(lon, lat, latlon=True, s=1, alpha=0.5, zorder=2)


        plt.title('Estimated number of clusters: %d' % n_clusters)

        #Upload an image to Cloud Storage
        plt.savefig(constants.C_MAP_FPATH)

        if os.path.exists(constants.C_MAP_FPATH):
            print "HERERE"
        else:
            print "NOT HEEEEEEEEEEEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRREEEEEEEEEEEE"

def setProjection(self,gid,axs=None,west=None,north=None,east=None,south=None):
        i = gid-1

        #If there is input on grid extent change the corners
        if (west != None and north != None and east != None and south != None):
            llcrnrlon = west
            llcrnrlat = south
            urcrnrlon = east
            urcrnrlat = north
        else:
            llcrnrlon = self.ll_lon[i]
            llcrnrlat = self.ll_lat[i]
            urcrnrlon = self.ur_lon[i]
            urcrnrlat = self.ur_lat[i]        


    #GOES is not in a specific projection so using cyl    
#        self.map[i] = Basemap(ax=axs,projection=self.projectionType,
#                                   llcrnrlon=self.ll_lon[i],llcrnrlat=self.ll_lat[i],
#                                   urcrnrlat = self.ur_lat[i],urcrnrlon = self.ur_lon[i],
#                                   resolution=self.resolution)
        self.map[i] = Basemap(ax=axs, projection=self.projectionType,
                      lat_ts = self.lat0[i], llcrnrlon = llcrnrlon,
                      llcrnrlat = llcrnrlat, urcrnrlon = urcrnrlon,
                      urcrnrlat = llcrnrlon, resolution=self.resolution)

        self.maplon[i],self.maplat[i] = self.map[i].makegrid(self.nx[i],self.ny[i])

def setProjection(self,gid=None,axs=None,west=None,north=None,east=None,south=None):
        self.map = [None]*1
        lon = self.variables['lon']
        lat = self.variables['lat']
        self.glons[0], self.glats[0] = np.meshgrid(lon,lat)
        self.nx = [self.glons[0].shape[1]]
        self.ny = [self.glons[0].shape[0]]
        self.dx = [50000.]
        self.dy = [62500.]
        if self.projectionType == None or self.projectionType == 'robin':
            self.projectionType = 'robin'
            self.lon0 = [0.0]
            self.lat0 = [0.0]
            self.ll_lon = [-179.375]
            self.ll_lat = [-90]
            self.ur_lon = [180]
            self.ur_lat = [90]
            self.map[0] = Basemap(ax=axs, projection=self.projectionType, 
                          lat_0 = self.lat0[0],
                      lon_0 = self.lon0[0], llcrnrlon = self.ll_lon[0],
                      llcrnrlat = self.ll_lat[0], urcrnrlon = self.ur_lon[0],
                      urcrnrlat = self.ur_lat[0], resolution=self.resolution)
        else:
            self.lon0 = [270.0]
            self.lat0 = [0.0]
            self.map[0] = Basemap(ax=axs, projection=self.projectionType,
                          lon_0 = self.lon0[0], boundinglat = self.lat0[0],
                          resolution=self.resolution)

def setProjection(self,gid,axs=None,west=None,north=None,east=None,south=None):
        i = gid-1

        #If there is input on grid extent change the corners
        if (west != None and north != None and east != None and south != None):
            llcrnrlon = west
            llcrnrlat = south
            urcrnrlon = east
            urcrnrlat = north
        else:
            llcrnrlon = self.ll_lon[i]
            llcrnrlat = self.ll_lat[i]
            urcrnrlon = self.ur_lon[i]
            urcrnrlat = self.ur_lat[i]

    #GOES is not in a specific projection so using cyl    
#        self.map[i] = Basemap(ax=axs,projection=self.projectionType,
#                                   llcrnrlon=self.ll_lon[i],llcrnrlat=self.ll_lat[i],
#                                   urcrnrlat = self.ur_lat[i],urcrnrlon = self.ur_lon[i],
#                                   resolution=self.resolution)
        self.map[i] = Basemap(ax=axs, projection=self.projectionType,
                      lat_ts = self.lat0[i], llcrnrlon = llcrnrlon,
                      llcrnrlat = llcrnrlat, urcrnrlon = urcrnrlon,
                      urcrnrlat = urcrnrlat, resolution=self.resolution)

        self.maplon[i],self.maplat[i] = self.map[i].makegrid(self.nx[i],self.ny[i])

    #This function navigate the input GOES file
    # i.e. determines the lat/long for each pixel
    # Based on the GOES-R Product Definition and User's Guide (PUG)

def setProjection(self,gid=None,axs=None,west=None,north=None,east=None,south=None):

        #If there is input on grid extent change the corners
        if (west != None and north != None and east != None and south != None):
            llcrnrlon = west
            llcrnrlat = south
            urcrnrlon = east
            urcrnrlat = north
        else:
            #Hardcode for U.S. for now
            self.ll_lon = -122.0
            self.ll_lat = 20.0
            self.ur_lon = -63.0
            self.ur_lat = 50.0
            llcrnrlon = self.ll_lon
            llcrnrlat = self.ll_lat
            urcrnrlon = self.ur_lon
            urcrnrlat = self.ur_lat

        self.map = [None]*1
        self.nx = [None]
        self.ny = [None]
        #Hardcode for U.S.
        self.lat0 = [38.0]
        self.lon0 = [-95.0] 
        self.map[0] = Basemap(projection=self.projectionType,
                      lat_0 = self.lat0[0],
                      lon_0 = self.lon0[0], llcrnrlon = llcrnrlon,
                      llcrnrlat = llcrnrlat, urcrnrlon = urcrnrlon,
                      urcrnrlat = urcrnrlat, resolution=self.resolution,
                      ax = axs)

############# GUI Related Functions ##################

    #Change the year

    #Sets the year

def setProjection(self,gid=None,axs=None,west=None,north=None,east=None,south=None):
        self.map = [None]*1
        self.nx = [None]
        self.ny = [None]
        self.glons = [None]*1
        self.glats = [None]*1
        self.ll_lon = self.lon0[0] - 2.25
        self.ll_lat = self.lat0[0] - 1.75
        self.ur_lon = self.lon0[0] + 2.25
        self.ur_lat = self.lat0[0] + 1.75
        self.map[0] = Basemap(projection=self.projectionType,
                      lat_0 = self.lat0[0],
                      lon_0 = self.lon0[0], llcrnrlon = self.ll_lon,
                      llcrnrlat = self.ll_lat, urcrnrlon = self.ur_lon,
                      urcrnrlat = self.ur_lat, resolution=self.resolution,
                      ax = axs)

        display = pyart.graph.RadarMapDisplay(self.radar)
        if ( self.var == 'velocity' or self.var == 'spectrum_width' ) and self.currentSweep == 0:
            self.currentSweep += 1
        if (self.var == 'differential_reflectivity' or 
            self.var == 'differential_phase' or 
            self.var == 'cross_correlation_ratio') and self.currentSweep == 1:
            self.currentSweep -= 1
        x,y = display._get_x_y(self.currentSweep,True,None)
        x0,y0 = self.map[0](self.lon0[0],self.lat0[0])
        self.glons[0],self.glats[0] = self.map[0]((x0+x*1000.),(y0+y*1000.),inverse=True)

def Mapa(self):
        self.mapa =  Basemap(projection='mill',lat_ts=10,llcrnrlon=self.lon0, urcrnrlon=self.lon1, llcrnrlat=self.lat0, urcrnrlat=self.lat1, resolution=self.res)

def build_map_obj(shapefile):
    shp = fiona.open(shapefile + '.shp')
    bds = shp.bounds

    shp.close()
    extra = 0.01
    ll = (bds[0], bds[1])
    ur = (bds[2], bds[3])
    coords = list(chain(ll, ur))

    w, h = coords[2] - coords[0], coords[3] - coords[1]

    m = Basemap(
                projection='tmerc',
                lon_0=(coords[0] + coords[2])/2,
                lat_0=(coords[1] + coords[3])/2,
                ellps='helmert',
                llcrnrlon=coords[0] - w * 0.01,
                llcrnrlat=coords[1] - extra + 0.01 * h,
                urcrnrlon=coords[2] + w * 0.01,
                urcrnrlat=coords[3] + extra + 0.01 * h,
                resolution='i',
                suppress_ticks=True
                )

    m.readshapefile(shapefile,
                    'SF',
                    color='black',
                    zorder=2
                    )
    return m, coords

def basemap(self, geobounds, **kwargs):
        """Return a :class:`matplotlib.mpl_toolkits.basemap.Basemap` object 
        for the map projection.

        Args:

            geobounds (:class:`wrf.GeoBounds`, optional):  The geobounds to 
                get the extents.  If set to None and using the *var* parameter,
                the geobounds will be taken from the variable.  If using a 
                file, then the geobounds will be taken from the native grid.

            **kwargs: Keyword arguments for creating a 
                :class:`matplotlib.mpl_toolkits.basemap.Basemap`.  By default,
                the domain bounds will be set to the native projection, the 
                resolution will be set to 'l', and the other projection 
                parameters will be set by the information in the file.

        Returns:

            :class:`matplotlib.mpl_toolkits.basemap.Basemap`: A Basemap
            object for the projection.

        See Also:

            :class:`matplotlib.mpl_toolkits.basemap.Basemap`

        """
        if not basemap_enabled():
            raise RuntimeError("'mpl_toolkits.basemap' is not "
                               "installed or is disabled")

        return self._basemap(geobounds, **kwargs)

def _proj4(self):
        _proj4 = ("+proj=eqc +units=meters +a={} +b={} "
                       "+lon_0={}".format(Constants.WRF_EARTH_RADIUS,
                                          Constants.WRF_EARTH_RADIUS,
                                          self.stand_lon))
        return _proj4

# Notes (may not be correct since this projection confuses me):
# Each projection system handles this differently.
# 1) In WRF, if following the WPS instructions, POLE_LON is mainly used to 
#    determine north or south hemisphere.  In other words, it determines if 
#    the globe is tipped toward or away from you.
# 2) In WRF, POLE_LAT is always positive, but should be negative in the 
#    proj4 based systems when using the southern hemisphere projections.
# 3) In cartopy, pole_longitude is used to describe the dateline, which 
#    is 180 degrees away from the normal central (standard) longitude 
#    (e.g. center of the projection), according to the cartopy developer.  
# 4) In basemap, lon_0 should be set to the central (standard) longitude.
# 5) In either cartopy, basemap or pyngl, I'm not sure that projections with
#    a pole_lon not equal to 0 or 180 can be plotted.  Hopefully people 
#    follow the WPS instructions, otherwise I need to see a sample file.
# 6) For items in 3 - 4, the "longitude" (lon_0 or pole_longitude) is 
#    determined by WRF's 
#    STAND_LON values, with the following calculations based on hemisphere:
#    BASEMAP:  NH:  -STAND_LON;  SH:  180.0 - STAND_LON
#    CARTOPY:  NH:  -STAND_LON - 180.; SH:  -STAND_LON
# 9) For PYNGL/NCL, you only need to set the center lat and center lon,
#    Center lat is the offset of the pole from +/- 90 degrees.  Center 
#    lon is -STAND_LON in NH and 180.0 - STAND_LON in SH.  
# 10) It also appears that NetCDF CF has no clear documentation on what 
#    each parameter means.  Going to assume it is the same as basemap, since
#    basemap appears to mirror the WMO way of doing things (tilt earth, then 
#    spin globe).
# 11) Basemap and cartopy produce projections that differ in their extent 
#     calculations by either using negative values or 0-360 (basemap).  For
#     this reason, the proj4 string for this class will use cartopy's values
#     to keep things in the -180 to 180, -90 to 90 range.
# 12) This projection makes me sad.

def __init__(self, *args, **kwargs):
        """Initialisation of the map object

        :param *args: additional non-keyword parameters (passed to `basemap 
            <http://matplotlib.org/basemap/api/basemap_api.html#mpl
            _toolkits.basemap.Basemap>`_)
        :param **kwargs: additional keyword parameters (passed to `basemap 
            <http://matplotlib.org/basemap/api/basemap_api.html#mpl
            _toolkits.basemap.Basemap>`_)

        .. note::

            Additional possible input in **kwargs wit to ``Basemap`` 
            objects is "topo_data" which will be set in case it is a valid 
            input (i.e. :class:`TopoData`) and will be used for plotting 
            topography
        """
        super(Map, self).__init__(*args, **kwargs)

        self.topo_data = None
        #IMPORTANT HANDLES
        self.contour_lines = None
        self.contour_filled = None

        self.colorbars = {}
        self.points = {}
        self.lines = {}
        self.polygons = {}
        self.texts = {}
        self.plotted_data = {}

        self.map_scale = None
        self.meridians = None
        self.parallels = None

        self.default_colors = {"contour_lines"  :   "#708090",
                               "land"           :   "#FFE0B2",
                               "water"          :   "#e6e6fa"}

def __init__(self, ):
        self.m = Basemap()

def map_world(self, res='l'):
        """
        It creates a map of the world.
        :return:
        """
        self.m = Basemap(projection='mill', resolution=res)
        # Draw the costal lines
        self.m.drawcoastlines()
        # Fill the continents with black
        self.m.fillcontinents(color='K')

def map_mediterranean(self, res='l'):
        """
        It creates a map of the Mediterranean.
        :param res: Resolution of the map
        :return:
        """
        self.m = Basemap(projection='mill', resolution=res, llcrnrlat=30, llcrnrlon=-13, urcrnrlat=47, urcrnrlon=38)
        # Draw the costal lines
        self.m.drawcoastlines()
        # Fill the continents with black
        self.m.fillcontinents(color='K')

def map(dataframe, title = "Map", colorbarName = None):

    fig, ax = plt.subplots(figsize=(80,40))

    m = Basemap(resolution='l', # c, l, i, h, f or None
            projection='robin',
            lon_0=0)

    m.drawmapboundary(fill_color='#46bcec')
    m.fillcontinents(color='#f2f2f2',lake_color='#46bcec')
    # m.drawcoastlines()

    plt.title(title, fontsize=50, y=1.08)

    m.readshapefile('visualization/World/World', 'world',drawbounds=False)

    df_plot = pd.DataFrame({
        'shapes': [Polygon(np.array(shape), True) for shape in m.world],
        'country': [country['ISO3'] for country in m.world_info]
        })

    df_plot = df_plot.merge(dataframe, on='country', how='left')

    df_plot = df_plot.dropna()

    cmap = plt.get_cmap('RdYlGn')   
    pc = PatchCollection(df_plot.shapes, zorder=2)
    norm = Normalize()

    pc.set_facecolor(cmap(norm(df_plot['value'].values)))
    ax.add_collection(pc)

    mapper = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)

    mapper.set_array(df_plot['value'])
    cbar = plt.colorbar(mapper, shrink=0.7, label = colorbarName)


    fig = plt.gcf()
    fig.savefig("Map.jpg")

    plt.show()

def draw_training_points(self, filename='points.pdf', world=False, figsize=(4,3)):
        '''
        draws training points on map
        '''
        import matplotlib as mpl
        mpl.use('Agg')
        import matplotlib.patches as mpatches
        import matplotlib.pyplot as plt
        from mpl_toolkits.basemap import Basemap, cm, maskoceans


        fig = plt.figure(figsize=figsize)
        lllat = 24.396308
        lllon = -124.848974
        urlat =  49.384358
        urlon = -66.885444
        if world:
            lllat = -90
            lllon = -180
            urlat = 90
            urlon = 180
        m = Basemap(llcrnrlat=lllat,
        urcrnrlat=urlat,
        llcrnrlon=lllon,
        urcrnrlon=urlon,
        resolution='c', projection='cyl')
        m.drawmapboundary(fill_color = 'white')
        m.drawcoastlines(linewidth=0.2)
        m.drawcountries(linewidth=0.2)

        ax = plt.gca()
        ax.xaxis.set_visible(False) 
        ax.yaxis.set_visible(False) 
        for spine in ax.spines.itervalues(): 
            spine.set_visible(False) 

        #fig = plt.figure()  # figsize=(4,4.2)

        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.spines['bottom'].set_visible(False)
        ax.spines['left'].set_visible(False)
        train_locs = self.df_train[['lat', 'lon']].values
        mlon, mlat = m(*(train_locs[:,1], train_locs[:,0]))
        #m.scatter(mlon, mlat, color='red', s=0.6)
        m.plot(mlon, mlat, 'r.', markersize=1)
        m.drawlsmask(land_color='lightgray',ocean_color="#b0c4de", lakes=True)
        plt.tight_layout()
        plt.savefig(filename)
        plt.close()
        print("the plot saved in " + filename)

def plot_map(outdict, e, savename, fig1=None, m=None):
    """
        This function will plot the output data in a scatter plot over a map of the
        satellite path.

        Args:
            outdict (dict[str, obj]): Output dictionary from analyzebeacons.
            e (dict[str, obj]): Output dictionary from ephem_doponly.
            savename(:obj:`str`): Name of the file the image will be saved to.
            fig1(:obj:'matplotlib figure'): Figure.
            m (:obj:'basemap obj'): Basemap object
    """

    t = outdict['time']
    slat = e['site_latitude']
    slon = e['site_longitude']
    if fig1 is None:
        fig1 = plt.figure()
    plt.figure(fig1.number)
    latlim = [math.floor(slat-15.), math.ceil(slat+15.)]
    lonlim = [math.floor(slon-15.), math.ceil(slon+15.)]
    if m is None:
        m = Basemap(lat_0=slat, lon_0=slon, llcrnrlon=lonlim[0], llcrnrlat=latlim[0],
                    urcrnrlon=lonlim[1], urcrnrlat=latlim[1])

    m.drawcoastlines(color="gray")

    m.plot(slon, slat, "rx")
    scat = m.scatter(e["sublon"](t), e["sublat"](t), c=outdict['rTEC'], cmap='viridis',
              vmin=0, vmax=math.ceil(sp.nanmax(outdict['rTEC'])))
    plt.title('Map of TEC Over Satellite Path')
    cb = plt.colorbar(scat,label='rTEC in TECu')
    fig1.savefig(savename, dpi=300)



    #plt.draw()

    scat = m.scatter(e["sublon"](t), e["sublat"](t), c=outdict['rTEC_amp'], cmap='viridis',
              vmin=0, vmax=math.ceil(sp.nanmax(outdict['rTEC_amp'])))
    plt.title('Map of TEC_Amp Over Satellite Path')
    cb.set_clim(vmin=0, vmax=math.ceil(sp.nanmax(outdict['rTEC_amp'])))
    cb.draw_all()

    #plt.tight_layout()
    figpath,name = os.path.split(savename)
    savename = os.path.join(figpath,'amp'+name)
    fig1.savefig(savename, dpi=300)
    plt.close(fig1)

#%% I/O for measurements

def Frame(self, grib_object, file_name, vmin, vmax):
        """ Initialize the first image of GIF
        grib_object: an object containing raw data to be visualized
        vmin: min of all data
        vmax: max of all data
        return generated figure instance
        """
        fig = plt.figure(figsize=(8,8))
        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

        lat,lon = grib_object.latlons()
        unit = grib_object['units']
        data_type = grib_object['name']
        date = self.GetTime(grib_object)
        data = grib_object.values

        m = Basemap(
                resolution='c', # c, l, i, h, f or None
                projection='cyl',
                lat_0=39.72, lon_0=-86.29,
                llcrnrlon=-87.79, llcrnrlat= 38.22,
                urcrnrlon=-84.79, urcrnrlat=41.22)

        parallels = np.arange(38.22, 41.22, 0.5)
        m.drawparallels(parallels,labels=[False,True,True,False])
        meridians = np.arange(-87.79, -84.79, 0.5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        x,y = m(lon, lat)
        im = m.pcolormesh(x,y,data,
                        shading='flat',
                        vmin=vmin,
                        vmax=vmax,
                        cmap=plt.cm.jet)
        # self.im = m.imshow(grib_object['values'],vmin=vmin,vmax=vmax, cmap=plt.cm.jet)
        cbar = plt.colorbar(location='bottom', fraction=0.046, pad=0.06)

         # Adjust the position of Unit
        cbar_ax = cbar.ax
        cbar_ax.text(0.0, -1.3, unit, horizontalalignment='left')
        m.readshapefile(self.shapeFile,'areas')
        title = plt.title(data_type + '   ' + date, fontsize = 'x-large')
        plt.savefig(file_name)
        plt.close(fig)

def plot_pierce_points(self,depth=410.0,ax='None',**kwargs):
   ####################################################################################
      '''
      Plots pierce points for a given depth.  If an axes object is supplied
      as an argument it will use it for the plot.  Otherwise, a new axes
      object is created.

      kwargs:
             depth:         pierce point depth
             proj:          map projection.  if none given, a default basemap axis is
                            made, defined by the coordinates of the corners. if 'ortho'
                            is given, a globe centered on Lat_0 and Lon_0 is made.
             ll_corner_lat: latitude of lower left corner of map
             ll_corner_lon: longitude of lower left corner of map
             ur_corner_lat: latitude of upper right corner of map
             Lat_0:         latitude center of ortho
             Lon_0:         logitude center of ortho
             return_ax:     whether or you return the basemap axis object. default False

      '''
      proj=kwargs.get('proj','default')
      ll_lat=kwargs.get('ll_corner_lat',-35.0)
      ll_lon=kwargs.get('ll_corner_lon',0.0)
      ur_lat=kwargs.get('ur_corner_lat',35.0)
      ur_lon=kwargs.get('ur_corner_lon',120.0)
      Lat_0=kwargs.get('Lat_0',0.0)
      Lon_0=kwargs.get('Lon_0',0.0)
      return_ax =kwargs.get('return_ax',False)

      if ax == 'None':
         if proj == 'default':
            m = Basemap(llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat)
         elif proj == 'ortho':
            m = Basemap(projection='ortho',lat_0=Lat_0,lon_0=Lon_0)
         m.drawmapboundary()
         m.drawcoastlines()
         m.fillcontinents()
      else:
         m = ax

      found_points = False
      for ii in self.pierce:
         if ii['depth'] == depth:
            x,y = m(ii['lon'],ii['lat'])
            found_points = True

      if found_points == True:
         m.scatter(x,y,100,marker='+',color='r',zorder=99)
      else:
         print "no pierce points found for the given depth"


   ###############################################################################

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

#########################################################################
# Plot multiple traces ontop of eachother, aligned on a phase and windowed
#########################################################################