Python mpl_toolkits.basemap.Basemap() Examples

def polar_stere(lon_w, lon_e, lat_s, lat_n, **kwargs):
        '''Returns a Basemap object (NPS/SPS) focused in a region.
        
        lon_w, lon_e, lat_s, lat_n -- Graphic limits in geographical coordinates.
                                      W and S directions are negative.
        **kwargs -- Aditional arguments for Basemap object.
        
        '''
        lon_0 = lon_w + (lon_e - lon_w) / 2.
        ref = lat_s if abs(lat_s) > abs(lat_n) else lat_n
        lat_0 = math.copysign(90., ref)
        proj = 'npstere' if lat_0 > 0 else 'spstere'
        prj = Basemap(projection=proj, lon_0=lon_0, lat_0=lat_0,
                              boundinglat=0, resolution='l')
        #prj = pyproj.Proj(proj='stere', lon_0=lon_0, lat_0=lat_0)
        lons = [lon_w, lon_e, lon_w, lon_e, lon_0, lon_0]
        lats = [lat_s, lat_s, lat_n, lat_n, lat_s, lat_n]
        x, y = prj(lons, lats)
        ll_lon, ll_lat = prj(min(x), min(y), inverse=True)
        ur_lon, ur_lat = prj(max(x), max(y), inverse=True)
        return Basemap(projection='stere', lat_0=lat_0, lon_0=lon_0,
                           llcrnrlon=ll_lon, llcrnrlat=ll_lat,
                           urcrnrlon=ur_lon, urcrnrlat=ur_lat, round=True,
                           resolution='l') 

def plot(self, variable, time, level, **plot_basemap_kwargs):
        """
        Wrapper to plot a specified field from an CFSReanalysis object.

        :param variable: str: variable to retrieve
        :param time: datetime: requested time
        :param level: int: requested pressure level
        :param plot_basemap_kwargs: kwargs passed to the plot.plot_functions.plot_basemap function (see the doc for
            plot_basemap for more information on options for Basemap plot)
        :return: matplotlib Figure object
        """
        from ..plot import plot_basemap
        print('CFSReanalysis.plot: plot of %s at %d mb (%s)' % (variable, level, time))
        field = self.field(variable, time, level)
        fig = plot_basemap(self.basemap, self.lon, self.lat, field, **plot_basemap_kwargs)
        return fig


# ==================================================================================================================== #
# CFSReanalysis object class
# ==================================================================================================================== # 

def _basemap_get_quicklook(area_def, data, vmin=None, vmax=None,
                           label='Variable (units)', num_meridians=45,
                           num_parallels=10, coast_res='110m', cmap='RdBu_r'):
    """Doing quicklook image plots with Basemap."""
    if area_def.shape != data.shape:
        raise ValueError('area_def shape %s does not match data shape %s' %
                         (list(area_def.shape), list(data.shape)))
    import matplotlib.pyplot as plt
    bmap = area_def2basemap(area_def, resolution=coast_res)
    bmap.drawcoastlines()
    if num_meridians > 0:
        bmap.drawmeridians(np.arange(-180, 180, num_meridians))
    if num_parallels > 0:
        bmap.drawparallels(np.arange(-90, 90, num_parallels))
    if not (np.ma.isMaskedArray(data) and data.mask.all()):
        col = bmap.imshow(data, origin='upper', vmin=vmin, vmax=vmax, cmap=cmap)
        plt.colorbar(col, shrink=0.5, pad=0.05).set_label(label)
    return plt 

def draw_geolines(area, dlon, dlat, basemap, linewidth=1):
    """
    Draw the parallels and meridians on a basemap plot.

    Parameters:

    * area : list
        ``[west, east, south, north]``, i.e., the area where the lines will
        be plotted
    * dlon, dlat : float
        The spacing between the lines in the longitude and latitude directions,
        respectively (in decimal degrees)
    * basemap : mpl_toolkits.basemap.Basemap
        The basemap used for plotting (see :func:`~geoist.vis.giplt.basemap`)
    * linewidth : float
        The width of the lines

    """
    west, east, south, north = area
    basemap.drawmeridians(numpy.arange(west, east, dlon), labels=[0, 0, 0, 1],
                          linewidth=linewidth)
    basemap.drawparallels(numpy.arange(south, north, dlat),
                          labels=[1, 0, 0, 0], linewidth=linewidth) 

def plot_route(individual):
    m = Basemap(projection='lcc', resolution=None,
                width=5E6, height=5E6,
                lat_0=-15, lon_0=-56)

    plt.axis('off')
    plt.title("Shortest Route")

    for i in range(0, len(individual.genes)):
        x, y = m(individual.genes[i].lng, individual.genes[i].lat)

        plt.plot(x, y, 'ok', c='r', markersize=5)
        if i == len(individual.genes) - 1:
            x2, y2 = m(individual.genes[0].lng, individual.genes[0].lat)
        else:
            x2, y2 = m(individual.genes[i+1].lng, individual.genes[i+1].lat)

        plt.plot([x, x2], [y, y2], 'k-', c='r') 

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 generate_map(output, lats=[], lons=[], wesn=None):
    """
    Using Basemap and the matplotlib toolkit, this function generates a map and
    puts a red dot at the location of every IP addresses found in the list.
    The map is then saved in the file specified in `output`.
    """
    print("Generating map and saving it to {}".format(output))
    if wesn:
        wesn = [float(i) for i in wesn.split('/')]
        m = Basemap(projection='cyl', resolution='l',
                llcrnrlon=wesn[0], llcrnrlat=wesn[2],
                urcrnrlon=wesn[1], urcrnrlat=wesn[3])
    else:
        m = Basemap(projection='cyl', resolution='l')
    m.bluemarble()
    x, y = m(lons, lats)
    m.scatter(x, y, s=1, color='#ff0000', marker='o', alpha=0.3)
    plt.savefig(output, dpi=300, bbox_inches='tight') 

def draw_countries(basemap, linewidth=1, style='dashed'):
    """
    Draw the country borders using the given basemap.

    Parameters:

    * basemap : mpl_toolkits.basemap.Basemap
        The basemap used for plotting (see :func:`~geoist.vis.giplt.basemap`)
    * linewidth : float
        The width of the lines
    * style : str
        The style of the lines. Can be: 'solid', 'dashed', 'dashdot' or
        'dotted'

    """
    lines = basemap.drawcountries(linewidth=linewidth)
    lines.set_linestyles(style) 

def draw_coastlines(basemap, linewidth=1, style='solid'):
    """
    Draw the coastlines using the given basemap.

    Parameters:

    * basemap : mpl_toolkits.basemap.Basemap
        The basemap used for plotting (see :func:`~geoist.vis.giplt.basemap`)
    * linewidth : float
        The width of the lines
    * style : str
        The style of the lines. Can be: 'solid', 'dashed', 'dashdot' or
        'dotted'

    """
    lines = basemap.drawcoastlines(linewidth=linewidth)
    lines.set_linestyles(style) 

def initBasemap(self):
        # NOTE:  http://matplotlib.org/basemap/api/basemap_api.html
        # NOTE:  https://gist.github.com/dannguyen/eb1c4e70565d8cb82d63
        self.mapAxes.clear()

        basemap_kwargs = {}
        basemap_kwargs.update(self.DEFAULT_BASEMAP_KWARGS)
        basemap_kwargs.update(
            ax=self.mapAxes
        )
        self.basemap = Basemap(**basemap_kwargs)

        self.basemap.bluemarble(scale=0.25, alpha=0.42)
        self.basemap.drawcoastlines(color='#555566', linewidth=1)
        self.basemap.drawmeridians(np.arange(0, 360, 30))
        self.basemap.drawparallels(np.arange(-90, 90, 30))

        self.canvas.draw_idle() 

def get_map_format():
        """
        Source:
            https://matplotlib.org/basemap/api/basemap_api.html
        Returns:
            Basemap: Constructed world basemap
        """
        world_map = Basemap(projection="cyl", llcrnrlat=-90, urcrnrlat=90,
                            llcrnrlon=-180, urcrnrlon=180, resolution="c")
        Plotter.draw_map_details(world_map)
        return world_map 

def BasePlot(self):

    plt.figure(figsize = (self._cfg['FigSize'][0],
                          self._cfg['FigSize'][1]))

    # Basemap
    self._map = Basemap(self._cfg['Bounds'][0],
                        self._cfg['Bounds'][1],
                        self._cfg['Bounds'][2],
                        self._cfg['Bounds'][3],
                        resolution = 'l',
                        projection = 'tmerc',
                        epsg = 3857)

    # Background land
    if self._cfg['Background'][0] == 'color':
      self._map.drawlsmask(land_color = self._cfg['Background'][1],
                           ocean_color = self._cfg['Background'][2],
                           grid = 1.25,
                           lakes = True)

    if self._cfg['Background'][0] == 'etopo':
      self._map.etopo(zorder = self._zo)

    if self._cfg['Background'][0] == 'esri':
      self._map.arcgisimage(service = self._cfg['Background'][1],
                            xpixels = self._cfg['Background'][2],
                            dpi = 300,
                            zorder = self._zo)

    if self._cfg['Background'][0] == 'relief':
      self._map.shadedrelief()

  #--------------------------------------------------------------------------------------- 

def map2DGrid(ax, grid, tstr, xlen=1.0, ylen=1.0, isLeft=False):
    """
    grid is a Grid2D object 
    """
    xmin,xmax,ymin,ymax = grid.getBounds()
    pdata = grid.getData()
    nr,nc = pdata.shape
    lonrange = np.linspace(xmin,xmax,num=nc)
    latrange = np.linspace(ymin,ymax,num=nr)
    lon,lat = np.meshgrid(lonrange,latrange)
    latmean = np.mean([ymin,ymax])
    lonmean = np.mean([xmin,xmax])
    m = Basemap(llcrnrlon=xmin,llcrnrlat=ymin,urcrnrlon=xmax,urcrnrlat=ymax,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='c',area_thresh=1000.,projection='lcc',\
            lat_1=latmean,lon_0=lonmean,ax=ax)
    # draw coastlines and political boundaries.
    m.drawcoastlines()
    #m.drawcountries()
    #m.drawstates()
    lons = np.arange(xmin,xmax,xlen)
    lats = np.arange(ymin,ymax,ylen)
    if isLeft:
        labels = labels=[1,0,0,0]
    else:
        labels = labels=[0,0,0,0]
    m.drawparallels(lats,labels=labels,color='white',fmt='%.1f') # draw parallels
    m.drawmeridians(lons,labels=[0,0,0,1],color='white',fmt='%.1f') # draw meridians
    pmesh = m.pcolormesh(lon,lat,np.flipud(grid.getData()),latlon=True)
    #plt.hold(True)
    ax.set_title(tstr)
    m.colorbar(pmesh) 

def plot_map(self):
        PC=np.loadtxt(DIR + '/hotspots.xy')
        PB=np.loadtxt(DIR + '/plate_boundaries.xy')


        self.fig2_cross = self.ui.mapfig2_cross.fig

        self.ax2_cross = self.fig2_cross.add_axes([0.01, 0.01, .98, .98])

        self.map = Basemap(projection='moll', lon_0=0, resolution="c",ax=self.ax2_cross)

        self.map.drawmapboundary(fill_color='#cccccc')
        self.map.fillcontinents(color='white', lake_color='#cccccc',zorder=0)

        self.fig2_cross.patch.set_alpha(0.0)

        self.fig2_cross.canvas.mpl_connect('button_press_event', self._on_map_mouse_click_event)
        # self.fig2_cross.canvas.mpl_connect('scroll_event',self._zoom_fun)

        xPC, yPC = self.map(PC[:,0],PC[:,1])
        xPB, yPB = self.map(PB[:,0],PB[:,1])
        self.__points_chauds = self.map.scatter(xPC, yPC, s=10, zorder=10,
                                                   color="magenta", marker="o",
                                                   edgecolor="k")
        colors=["7fff00"]
        self.__plate_boundaries = self.map.scatter(xPB, yPB, s=2, zorder=10,
                                                   color="0.3", marker=".")

        EQ=np.loadtxt(DIR + '/catalogue.xy')
        xEQ, yEQ = self.map(EQ[:,0],EQ[:,1])
        self.__eq_map_obj = self.map.scatter(xEQ, yEQ, s=10, zorder=10,
                                                       color="white", marker="o",
                                                       edgecolor="k")

        self.fig2_cross.canvas.draw() 

def test_area_def2basemap(self):
        """Test the area to Basemap object conversion function."""
        from pyresample import plot
        from pyresample import parse_area_file
        area_def = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml'), 'ease_sh')[0]
        bmap = plot.area_def2basemap(area_def)
        self.assertTrue(bmap.rmajor == bmap.rminor and bmap.rmajor == 6371228.0,
                        'Failed to create Basemap object') 

def show_quicklook(area_def, data, vmin=None, vmax=None,
                   label='Variable (units)', num_meridians=45,
                   num_parallels=10, coast_res='110m', cmap='RdBu_r'):
    """Display default quicklook plot.

    Parameters
    ---------
    area_def : object
        geometry.AreaDefinition object
    data : numpy array | numpy masked array
        2D array matching area_def. Use masked array for transparent values
    vmin : float, optional
        Min value for luminescence scaling
    vmax : float, optional
        Max value for luminescence scaling
    label : str, optional
        Label for data
    num_meridians : int, optional
        Number of meridians to plot on the globe
    num_parallels : int, optional
        Number of parallels to plot on the globe
    coast_res : {'c', 'l', 'i', 'h', 'f'}, optional
        Resolution of coastlines

    Returns
    -------
    bmap : Basemap object

    """
    plt = _get_quicklook(area_def, data, vmin=vmin, vmax=vmax,
                         label=label, num_meridians=num_meridians,
                         num_parallels=num_parallels, coast_res=coast_res,
                         cmap=cmap)
    plt.show()
    plt.close() 

def draw_bb_map(self, fig, ax, BB=None, range_rings=None,
                                lat_spacing=None, lon_spacing=None):
        '''Draw the Beam Blockage'''
        if BB is None:
            BB = self.CBB
            
        if lat_spacing is None:
            lat_spacing = 1.
        if lon_spacing is None:
            lon_spacing = 1.
            
        bm2 = Basemap(projection='cea', resolution='l', area_thresh = 10000.,
                    llcrnrlon=self.minlon, urcrnrlon=self.maxlon,
                    llcrnrlat=self.minlat, urcrnrlat=self.maxlat,
                    ax=ax)
        ax.set_title('Beam-blockage fraction', fontdict=TITLEDICT)
        bm2.drawmeridians(np.arange(self.minlon, self.maxlon, lon_spacing), labels=[1,0,0,1])           
        bm2.drawparallels(np.arange(self.minlat, self.maxlat, lat_spacing), labels=[1,0,0,1])             
        bm2.drawcountries()
        bm2.drawcoastlines()
        bm2.drawrivers()

        xbm2, ybm2 = bm2(self.rng_lon, self.rng_lat)
        BBmap = bm2.pcolormesh(xbm2, ybm2, BB, vmin=0., vmax=1., cmap=self.bb_cmap)
        if range_rings is not None:
            for nn in range(len(range_rings)):
                self.plot_range_ring(range_rings[nn], bm=bm2)
        fig.colorbar(BBmap, ax=ax) 

def plot_movie(m, lat, lon, val, pred, dates, model_title='', plot_kwargs=None, out_directory=None):
    """
    Plot a series of images for a forecast and the verification.

    :param m: Basemap object
    :param lat: ndarray (lat, lon): latitude values
    :param lon: ndarray (lat, lon): longitude values
    :param val: ndarray (t, lat, lon): verification
    :param pred: ndarray (t, lat, lon): predicted forecast
    :param dates: array-like: datetime objects of verification datetimes
    :param model_title: str: name of the model, e.g., 'Neural net prediction'
    :param plot_kwargs: dict: passed to the plot pcolormesh() method
    :param out_directory: str: folder in which to save image files
    """
    if (len(dates) != val.shape[0]) and (len(dates) != pred.shape[0]):
        raise ValueError("'val' and 'pred' must have the same first (time) dimension as 'dates'")
    plot_kwargs = plot_kwargs or {}
    fig = plt.figure()
    fig.set_size_inches(6, 4)
    x, y = m(lon, lat)
    dt = dates[1] - dates[0]
    for d, date in enumerate(dates):
        hours = (d + 1) * dt.total_seconds() / 60 / 60
        ax = plt.subplot(211)
        m.pcolormesh(x, y, val[d], **plot_kwargs)
        m.drawcoastlines()
        m.drawparallels(np.arange(0., 91., 45.))
        m.drawmeridians(np.arange(0., 361., 90.))
        ax.set_title('Verification (%s)' % date)
        ax = plt.subplot(212)
        m.pcolormesh(x, y, pred[d], **plot_kwargs)
        m.drawcoastlines()
        m.drawparallels(np.arange(0., 91., 45.))
        m.drawmeridians(np.arange(0., 361., 90.))
        ax.set_title('%s at $t=%d$ (%s)' % (model_title, hours, date))
        plt.savefig('%s/%05d.png' % (out_directory, d), bbox_inches='tight', dpi=150)
        fig.clear() 

def generate_basemap(self, llcrnrlat=None, llcrnrlon=None, urcrnrlat=None, urcrnrlon=None):
        """
        Generates a Basemap object for graphical plot of ERA5 data on a 2-D plane. Bounding box parameters
        are either given, or if None, read from the extremes of the loaded lat/lon data. Other projection parameters
        are set to the default ERA5 configuration.

        :param llcrnrlat: float: lower left corner latitude
        :param llcrnrlon: float: lower left corner longitude
        :param urcrnrlat: float: upper right corner latitude
        :param urcrnrlon: float: upper right corner longitude
        :return:
        """
        from mpl_toolkits.basemap import Basemap

        try:
            default = llcrnrlat * llcrnrlon * urcrnrlat * urcrnrlon  # error if any are None
            default = False
        except TypeError:
            default = True

        if default:
            try:
                lat = self.lat
                lon = self.lon
            except (AttributeError, KeyError):
                raise ValueError('I can generate a default Basemap with None parameters, but only if I have some '
                                 'data loaded first!')
            llcrnrlon, llcrnrlat = lon[0, 0], lat[-1, -1]
            urcrnrlon, urcrnrlat = lon[-1, -1], lat[0, 0]

        basemap = Basemap(projection='cyl', llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat,
                          llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon, resolution='l')

        self.basemap = basemap 

def draw_terrain_height_map(self, fig, ax, vmin=None, vmax=None,
                                lat_spacing=None, lon_spacing=None):
        '''Draw the terrain heights'''
        if vmin is None:
            topomin = 0.05
        else:
            topomin = vmin/1000.
        if vmax is None:
            topomax = 3.
        else:
            topomax = vmax/1000.
            
        if lat_spacing is None:
            lat_spacing = 1.
        if lon_spacing is None:
            lon_spacing = 1.
            
        bm1 = Basemap(projection='cea', resolution='l', area_thresh = 10000.,
                    llcrnrlon=self.minlon, urcrnrlon=self.maxlon,
                    llcrnrlat=self.minlat, urcrnrlat=self.maxlat,
                    ax=ax)
        ax.set_title('Terrain within %02d km of Radar (km)'%(self.range), fontdict=TITLEDICT)
        bm1.drawmeridians(np.arange(self.minlon, self.maxlon, lon_spacing), labels=[1,0,0,1])           
        bm1.drawparallels(np.arange(self.minlat, self.maxlat, lat_spacing), labels=[1,0,0,1])             
        bm1.drawcountries()
        bm1.drawcoastlines()
        bm1.drawrivers()
        
        xbm1, ybm1 = bm1(self.lon, self.lat)
        Htmap = bm1.pcolormesh(xbm1, ybm1, self.topo/1000., vmin=topomin, vmax=topomax,
                               cmap=self.terr_cmap)
        bm1.plot(self.rlon, self.rlat, 'rD', latlon=True)
        #plot_range_ring(50., bm=bm1, color='w')
        fig.colorbar(Htmap, ax=ax) 

def create_map(lon, lat):
    """
    Create a map projection.
    """
    lon_center = lon[0]
    lat_center = lat[0]
    return Basemap(
        lon_0=lon_center,
        lat_0=lat_center, projection='tmerc',
        width=1e-5, height=1e-5) 

def project_lat_lon(df):
    gps_map = Basemap(lat_0=df.GPS_Lat.values[0],
                      lon_0=df.GPS_Lon.values[0],
                      width=11e-5, height=1e-5, projection='tmerc')
    gps_y, gps_x = gps_map(df.GPS_Lon.values, df.GPS_Lat.values)
    gps_z = df.GPS_Alt - df.GPS_Alt.values[0]
    df_new = pandas.DataFrame(pandas.DataFrame({
        'GPS_X': gps_x, 'GPS_Y': gps_y, 'GPS_Z': gps_z}, index=df.index))
    return pandas.concat([df, df_new], axis=1)

# vim: set et fenc= ff=unix sts=0 sw=4 ts=4 : 

def basemap_from_file(ifile, withgrid=False, **kwds):
    """
    Typically, the user will need to provide some options
    """
    proj4 = getproj4(ifile, withgrid=withgrid)
    basemap_options = basemap_options_from_proj4(proj4, **kwds)
    if 'llcrnrx' in basemap_options:
        if 'urcrnrx' in kwds:
            basemap_options['urcrnrx'] = kwds['urcrnrx']
        elif 'width' in kwds:
            basemap_options['urcrnrx'] = basemap_options['llcrnrx'] + \
                kwds['width']
        elif 'x' in ifile.variables:
            x = ifile.variables['x']
            urx = x.max() + np.mean(np.diff(x))
            basemap_options['urcrnrx'] = urx
        else:
            raise KeyError(
                'When a false_easting is available, the file must contain ' +
                'an x variable or the user must supply width or urcrnrx')
    if 'llcrnry' in basemap_options:
        if 'urcrnry' in kwds:
            basemap_options['urcrnry'] = kwds['urcrnry']
        elif 'height' in kwds:
            basemap_options['urcrnry'] = basemap_options['llcrnry'] + \
                kwds['height']
        elif 'y' in ifile.variables:
            y = ifile.variables['y']
            ury = y.max() + np.mean(np.diff(y))
            basemap_options['urcrnry'] = ury
        else:
            raise KeyError(
                'When a false_northing is available, the file must contain ' +
                'a y variable or the user must supply height or urcrnry')

    from mpl_toolkits.basemap import Basemap
    print(basemap_options)
    bmap = Basemap(**basemap_options)
    return bmap 

def basemap_from_proj4(proj4, **kwds):
    from mpl_toolkits.basemap import Basemap
    basemap_options = basemap_options_from_proj4(proj4, **kwds)
    if basemap_options['projection'] in ('lonlat', 'longlat'):
        basemap_options['projection'] = 'cyl'
    bmap = Basemap(**basemap_options)
    return bmap 

def generate_basemap(self, llcrnrlat=None, llcrnrlon=None, urcrnrlat=None, urcrnrlon=None):
        """
        Generates a Basemap object for graphical plot of CFSR data on a 2-D plane. Bounding box parameters
        are either given, or if None, read from the extremes of the loaded lat/lon data. Other projection parameters
        are set to the default CFSR configuration.

        :param llcrnrlat: float: lower left corner latitude
        :param llcrnrlon: float: lower left corner longitude
        :param urcrnrlat: float: upper right corner latitude
        :param urcrnrlon: float: upper right corner longitude
        :return:
        """
        from mpl_toolkits.basemap import Basemap

        try:
            default = llcrnrlat * llcrnrlon * urcrnrlat * urcrnrlon  # error if any are None
            default = False
        except TypeError:
            default = True

        if default:
            try:
                lat = self.lat
                lon = self.lon
            except (AttributeError, KeyError):
                raise ValueError('I can generate a default Basemap with None parameters, but only if I have some '
                                 'data loaded first!')
            llcrnrlon, llcrnrlat = lon[0, 0], lat[-1, -1]
            urcrnrlon, urcrnrlat = lon[-1, -1], lat[0, 0]

        basemap = Basemap(projection='cyl', llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat,
                          llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon, resolution='l')

        self.basemap = basemap 

def imageMap2(lons, lats, vals, vmin=None, vmax=None, 
             imageWidth=None, imageHeight=None, outFile=None,
             projection='cyl', cmap=M.cm.jet, makeFigure=False,
             borders=[0., -90., 360., 90.], autoBorders=True, borderSlop=10.,
             meridians=[0, 360, 60], parallels=[-60, 90, 30],
             **options
             ):
#    lons = normalizeLons(lons)
    if vmin == 'auto': vmin = None
    if vmax == 'auto': vmax = None
    if imageWidth is not None: makeFigure = True
    if projection is None or projection == '': projection = 'cyl'
    if cmap is None or cmap == '': cmap = M.cm.jet
    if isinstance(cmap, types.StringType) and cmap != '':
        try:
            cmap = eval('M.cm.' + cmap)
        except:
            cmap = M.cm.jet

#    ensureItems(options, {'xlabel': 'Longitude (deg)', 'ylabel': 'Latitude (deg)', \
    ensureItems(options, { \
                     'title': 'An Image Map', 'dpi': 100,
                     'imageWidth': imageWidth or 1024, 'imageHeight': imageHeight or 768})
    if autoBorders:
        borders = [min(min(lons)), min(min(lats)), max(max(lons)), max(max(lats))]
        borders = roundBorders(borders, borderSlop)

    m = Basemap(borders[0], borders[1], borders[2], borders[3], \
                projection=projection, lon_0=N.average([lons.flat[0], lons.flat[-1]]))

    if makeFigure:
        dpi = float(options['dpi'])
        width = float(imageWidth) / dpi
        if imageHeight is None:
            height = width * m.aspect
        else:
            height = float(imageHeight) / dpi
        f = M.figure(figsize=(width,height)).add_axes([0.1,0.1,0.8,0.8], frameon=True)

    if vmin is not None or vmax is not None: 
        if vmin is None:
            vmin = min(min(vals))
        else:
            vmin = float(vmin)
        if vmax is None:
            vmax = max(max(vals))
        else:
            vmax = float(vmax)
        vrange = vmax - vmin
        levels = N.arange(vmin, vmax, vrange/30.)
    else:
        levels = 30

    c = m.contourf(lons, lats, vals, levels, cmap=cmap, colors=None)
    m.drawcoastlines()
    m.drawmeridians(range(meridians[0], meridians[1], meridians[2]), labels=[0,0,0,1])
    m.drawparallels(range(parallels[0], parallels[1], parallels[2]), labels=[1,1,1,1])
    M.colorbar(c, orientation='horizontal')
    evalKeywordCmds(options)
    if outFile: M.savefig(outFile, **validCmdOptions(options, 'savefig')) 

def stat_spatial(ifile0, ifile1, funcs=__all__,
                 variables=['O3'], counties=False):
    """
    ifile0 - left hand side of equation
    ifile1 - right hand side of equation
    variables - list of variables to plot
    """
    from mpl_toolkits.basemap import Basemap
    from matplotlib.pyplot import figure, show
    lonlatcoords = getattr(ifile0, 'lonlatcoords',
                           getattr(ifile1, 'lonlatcoords', ''))
    lon, lat = np.array(
        map(lambda x: map(float, x.split(',')), lonlatcoords.split('/'))).T
    latmin, latmax = lat.min(), lat.max()
    lonmin, lonmax = lon.min(), lon.max()
    bmap = Basemap(llcrnrlat=latmin, llcrnrlon=lonmin,
                   urcrnrlat=latmax, urcrnrlon=lonmax, resolution='i')
    for vark in variables:
        for statname in funcs:
            statfunc = eval(statname)
            fig = figure()
            ax = fig.add_subplot(111)
            ax.set_title(vark + ' ' + statname)
            var_0 = ifile0.variables[vark]
            var_1 = ifile1.variables[vark]
            try:
                pidx = list(var_0.dimensions).index('points')
            except Exception:
                pidx = list(var_1.dimensions).index('points')
            statvs = []
            for sitei in range(var_0.shape[pidx]):
                val_0 = var_0[:].take([sitei], axis=pidx).ravel()
                val_1 = var_1[:].take([sitei], axis=pidx).ravel()
                statvs.append(statfunc(val_0, val_1))
            dots = bmap.scatter(x=lon, y=lat, c=statvs, ax=ax, cmap='jet')

            bmap.drawcoastlines(ax=ax)
            bmap.drawcountries(ax=ax)
            bmap.drawstates(ax=ax)
            fig.colorbar(dots)
            if counties:
                bmap.counties(ax=ax)
            show() 

def render(d, lats, lons, z, primary, secondary, parameter):
    fig = plt.figure()
    ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

    ax.set_title(string.upper("%s vs. %s" % (primary, secondary)))
    # ax.set_ylabel('Latitude')
    # ax.set_xlabel('Longitude')

    minLatA = np.min(lats)
    maxLatA = np.max(lats)
    minLonA = np.min(lons)
    maxLonA = np.max(lons)

    minLat = minLatA - (abs(maxLatA - minLatA) * 0.1)
    maxLat = maxLatA + (abs(maxLatA - minLatA) * 0.1)

    minLon = minLonA - (abs(maxLonA - minLonA) * 0.1)
    maxLon = maxLonA + (abs(maxLonA - minLonA) * 0.1)

    minLon, maxLon, minLat, maxLat = __square(minLon, maxLon, minLat, maxLat)

    # m = Basemap(projection='mill', llcrnrlon=-180,llcrnrlat=-80,urcrnrlon=180,urcrnrlat=80,resolution='l')
    m = Basemap(projection='mill', llcrnrlon=minLon, llcrnrlat=minLat, urcrnrlon=maxLon, urcrnrlat=maxLat,
                resolution='l')

    m.drawparallels(np.arange(minLat, maxLat, (maxLat - minLat) / 5.0), labels=[1, 0, 0, 0], fontsize=10)
    m.drawmeridians(np.arange(minLon, maxLon, (maxLon - minLon) / 5.0), labels=[0, 0, 0, 1], fontsize=10)

    m.drawcoastlines()
    m.drawmapboundary(fill_color='#99ffff')
    m.fillcontinents(color='#cc9966', lake_color='#99ffff')

    # lats, lons = np.meshgrid(lats, lons)

    masked_array = np.ma.array(z, mask=np.isnan(z))
    z = masked_array

    values = np.zeros(len(z))
    for i in range(0, len(z)):
        values[i] = ((z[i] - np.min(z)) / (np.max(z) - np.min(z)) * 20.0) + 10

    x, y = m(lons, lats)

    im1 = m.scatter(x, y, values)

    im1.set_array(z)
    cb = m.colorbar(im1)

    units = PARAMETER_TO_UNITS[parameter] if parameter in PARAMETER_TO_UNITS else PARAMETER_TO_UNITS["sst"]
    cb.set_label("Difference %s" % units)

    sio = StringIO()
    plt.savefig(sio, format='png')
    plot = sio.getvalue()
    if d is not None:
        d['plot'] = plot
    return plot 

def MapPlot(
    data,
    fig,
    cmap=None,
    clim=None,
    bgcolor=np.array([1, 1, 1]) * 0.3,
    facecolor=np.array([1, 1, 1]) * 0.3,
    lons=None,
    lats=None,
    title=None,
    label=None,
    linewidth=None,
    norm=mpl.colors.Normalize(),
    resolution="c",
    proj="robin",
    lon_0=180,
):
    if lons is None:
        raise RuntimeError("map plotting needs lons input")
    if lats is None:
        raise RuntimeError("map plotting needs lats input")
    # cmap.set_bad(bgcolor, 1)
    # ax = fig.add_axes([0.15, 0.15, 0.85, 0.85])
    # ax.set_axis_off()
    # ax.set_facecolor(facecolor)
    # ax.set_aspect(1, anchor='C')
    m = Basemap(projection=proj, lon_0=lon_0, resolution=resolution)
    pixels = np.squeeze(np.ma.array(data, mask=np.isnan(data)))
    im = m.pcolor(
        np.array(lons),
        np.array(lats),
        pixels,
        cmap=cmap,
        vmin=clim[0],
        vmax=clim[1],
        norm=norm,
        linewidth=linewidth,
        latlon=True,
    )
    # TODO: Customize these eventually?
    m.drawmapboundary(fill_color=bgcolor, linewidth=1, color=bgcolor)
    m.drawcoastlines(color="0.75")
    m.fillcontinents(color="0.8")
    cb = m.colorbar(im, "right", size="3%", pad="8%")
    if label is not None:
        cb.set_label(label, fontsize=20, labelpad=10)
    cb.ax.tick_params(labelsize=20)
    if title is not None:
        plt.gca().set_title(title, y=1.03, fontsize=20) 

def tomomap_plot(model,para,depth,NSmax,Vmax,lon0=140,plate_boundaries=True,hotspots=True):

    path = DIR2 + "/output_files_map"
    dirs = os.listdir(path)
    os.chdir(DIR2 + '/output_files_map')
    for i in range(len(dirs)):
    	os.remove(dirs[i])
    os.chdir(DIR)
#    os.chdir('../../')

    Z2, lat, lon = tomomap(model,para,depth,NSmax)

    fig = plt.figure()
    model2 = Z2[:,2]
    X2, Y2 = np.meshgrid(lon, lat)
    model_new3 = np.reshape(model2, (len(lon),len(lat)))
    model_new4 = np.transpose(model_new3)

    ax  = fig.add_axes([0.1, 0.01, .8, .8])
    map = Basemap(projection='moll', lon_0=lon0, resolution="c",ax=ax)
    __s2 = map.pcolormesh(X2, Y2, model_new4, cmap='RdYlBu',latlon=True, vmin=-Vmax, vmax=Vmax)
    map.drawcoastlines(linewidth=1)
    step = float(Vmax)/4

    PC=np.loadtxt(DIR + '/hotspots.xy')
    PB=np.loadtxt(DIR + '/plate_boundaries.xy')
    xPC, yPC = map(PC[:,0],PC[:,1])
    xPB, yPB = map(PB[:,0],PB[:,1])
    if hotspots:
    	map.scatter(xPC, yPC, s=20, zorder=10,color="magenta", marker="o",edgecolor="k")
    if plate_boundaries:
    	map.scatter(xPB, yPB, s=2, zorder=10,color="gray", marker=".")

    path = DIR2 + "/output_files_map"
    dirs = os.listdir(path)
    os.chdir(DIR2 + '/output_files_map')
    for i in range(len(dirs)):
        os.remove(dirs[i])
    os.chdir(cwd)

    fig.colorbar(__s2,orientation="horizontal", ticks=np.arange(-Vmax,Vmax+step,step))
    plt.show()