def calc_ratio(fsoil_mary, fsoil_kettle):
    lon, lat, topo = sp.parse_STEM_coordinates(
        os.path.join(os.environ['SARIKA_INPUT'], 'TOPO-124x124.nc'))
    fsoil_mary = maskoceans(lon, lat, fsoil_mary)
    fsoil_kettle = maskoceans(lon, lat, fsoil_kettle)
    ratio = ma.masked_invalid(fsoil_kettle) / ma.masked_invalid(fsoil_mary)
    return(ratio)

def main(base_folder="/skynet3_rech1/huziy/veg_fractions/",
         fname="pm1983120100_00000000p", canopy_name="Y2C", label="USGS",
         depth_to_bedrock_name="8L"
         ):
    data_path = os.path.join(base_folder, fname)
    r = RPN(data_path)

    veg_fractions = r.get_2D_field_on_all_levels(name=canopy_name)
    print(list(veg_fractions.keys()))
    sand = r.get_first_record_for_name("SAND")
    clay = r.get_first_record_for_name("CLAY")

    dpth_to_bedrock = r.get_first_record_for_name(depth_to_bedrock_name)

    proj_params = r.get_proj_parameters_for_the_last_read_rec()

    lons, lats = r.get_longitudes_and_latitudes_for_the_last_read_rec()
    print(lons.shape)

    rll = RotatedLatLon(lon1=proj_params["lon1"], lat1=proj_params["lat1"],
                        lon2=proj_params["lon2"], lat2=proj_params["lat2"])

    lon0, lat0 = rll.get_true_pole_coords_in_rotated_system()
    plon, _ = rll.get_north_pole_coords()

    b = Basemap(projection="rotpole", llcrnrlon=lons[0, 0], llcrnrlat=lats[0, 0],
                urcrnrlon=lons[-1, -1], urcrnrlat=lats[-1, -1], lon_0=lon0 - 180,
                o_lon_p=lon0, o_lat_p=lat0)

    lons[lons > 180] -= 360
    for lev in list(veg_fractions.keys()):
        veg_fractions[lev] = maskoceans(lons, lats, veg_fractions[lev], inlands=False)

    sand = maskoceans(lons, lats, sand)
    clay = maskoceans(lons, lats, clay)
    dpth_to_bedrock = maskoceans(lons, lats, dpth_to_bedrock)

    x, y = b(lons, lats)
    plot_veg_fractions(x, y, b, veg_fractions, out_image=os.path.join(base_folder,
                                                                      "veg_fractions_{0}.jpeg".format(label)))
    plot_sand_and_clay(x, y, b, sand, clay, out_image=os.path.join(base_folder,
                                                                   "sand_clay_{0}.jpeg".format(label)))

    # set relation between vegetation frsction fields and names
    veg_fract_dict = {}
    for lev, the_field in veg_fractions.items():
        lev = int(lev)
        if lev not in y2c_level_to_title:
            continue
        veg_fract_dict[y2c_level_to_title[lev]] = the_field

    data = {
        "SAND": sand, "CLAY": clay, "BDRCK_DEPTH": dpth_to_bedrock
    }
    data.update(veg_fract_dict)

    return b, lons, lats, data, label

def get_mask(arr, lats=None, lons=None, masknan=True, maskocean=False, maskland=False):
    '''
        Return array which is a mask for the input array
        to mask the ocean or land you need to put in the lats, lons of the data
    '''
    mask=np.zeros(np.shape(arr),dtype=np.bool)
    if masknan:
        mask=np.isnan(arr)
    if maskocean or maskland:
        if len(np.shape(lats)) == 1:
            lons,lats = np.meshgrid(lons,lats)
    if maskocean:
        mask = mask + maskoceans(lons,lats,arr, inlands=False).mask
    if maskland:
        mask = mask + ~(maskoceans(lons,lats,arr, inlands=False).mask)
    return mask

def ResearchRegion_surface():
    """
    在地图上画出柱表面混合比图
    :return:
    """

    fig = plt.figure(figsize=(11, 8), facecolor="white")
    # data = np.loadtxt('seasonAvr_data/SurfaceMixingRatio/1_seasonAvr.txt')
    data = np.loadtxt("allYearAvr_data/SurfaceMixingRatio/allYearAvr.txt")
    arr = np.zeros((180, 360))
    for i in range(180):
        arr[i, :] = data[179 - i, :]

    longitude = np.loadtxt("lonlat_data/longitude.txt")
    latitude = np.loadtxt("lonlat_data/latitude.txt")

    m = Basemap(llcrnrlon=70, llcrnrlat=15, urcrnrlon=138, urcrnrlat=55, projection="mill", resolution="h")
    m.drawparallels(np.arange(5.5, 90.5, 1.0), color="w", linewidth=0.5, dashes=[1, 1], labels=[0, 0, 0, 0])
    m.drawmeridians(np.arange(60.5, 181.5, 1.0), color="w", linewidth=0.5, dashes=[1, 1], labels=[0, 0, 0, 0])
    m.drawmapboundary(fill_color="0.3")
    m.readshapefile("shp/CHINA", "CHINA", drawbounds=1, color="black")

    topo = maskoceans(longitude, latitude, arr)
    im = m.pcolormesh(longitude, latitude, topo, shading="flat", cmap=plt.cm.jet, latlon=True, vmin=0, vmax=500)
    m.drawlsmask(ocean_color="w", lsmask=0)

    cbar = m.colorbar()
    cbar.ax.set_ylabel("SurfaceMixingRatio", color="black", fontsize="14", rotation=90)
    plt.show()

def main():

    # create the image folder if necessary
    img_folder = "bulk_field_capacity_model"
    if not os.path.isdir(img_folder):
        os.mkdir(img_folder)

    data, lons, lats, bmp = get_data_and_coords()
    lons[lons > 180] -= 360
    print(list(data.keys()))

    # reproject coords
    x, y = bmp(lons, lats)

    clevs = np.arange(0, 0.5, 0.02)
    cmap = cm.get_cmap("rainbow", lut=len(clevs))
    # plot for all levels right away
    for lev, field in data.items():
        fig = plt.figure()
        plt.title(r"$\theta_{\rm fc}$, " + "soil lev = {}".format(lev))
        to_plot = maskoceans(lons, lats, field, inlands=True)
        img = bmp.contourf(x, y, to_plot, levels=clevs, cmap=cmap)
        bmp.colorbar(img)

        bmp.drawcoastlines()

        print("lev={}, fc-min={}, fc-max={}".format(lev, field.min(), field.max()))

        fname = "thfc_lev_{}.png".format(lev)
        fig.savefig(os.path.join(img_folder, fname))
        plt.close(fig)

def plot_for_simulation(axis=None, sim_path="", cmap=None, cnorm=None,
                        start_year=None, end_year=None, months=None):
    """
    plot a panel for each simulation

    :param axis:
    :param sim_path:
    :param cmap:
    :param cnorm:
    """

    if months is None:
        months = list(range(1, 13))

    lons, lats, bm = analysis.get_basemap_from_hdf(sim_path)

    params = dict(
        path1=sim_path, path2=sim_path,
        start_year=start_year, end_year=end_year,
        varname1="I1", level1=0,
        varname2="AV", level2=0,
        months=months
    )

    corr, i1_clim, av_clim = calculate_correlation_field_for_climatology(**params)

    # convert longitudes to the [-180, 180] range
    lons[lons > 180] -= 360
    corr = maskoceans(lons, lats, corr)

    x, y = bm(lons, lats)

    im = bm.pcolormesh(x, y, corr, norm=cnorm, cmap=cmap, ax=axis)
    bm.drawcoastlines()
    return im, corr

def latlonggrid(longmin, longmax, latmin, latmax, step, **kwargs):
    """ Produce a dense grid suitable for drawing climate maps -- land areas only.
    
    The number of grid points per degree is the same for longitude and latitude at the equator.
    Away from the equator, we thin it out a bit more to keep the number of grid points
    per kilometer roughhly constant.
    
    Parameters --
    longmin -- minimum longtiude
    longmax -- maximum longtiude
    latmin -- minimum latitude
    latmax -- maximum latitude
    step -- The number of grid points per degree for latitude and for longitude at the equator.
    """
    llgrid = pd.DataFrame(
        [
            [lat / step, long / step / cos(pi * lat / step / 180)]
            for lat in range(latmin * step, latmax * step + 1)
            for long in range(
                floor(longmin * step * cos(pi * lat / step / 180)), floor(longmax * step * cos(pi * lat / step / 180))
            )
        ],
        columns=["lat", "long"],
    )

    masked = maskoceans(llgrid.long, llgrid.lat, llgrid.index.get_level_values(0), **kwargs)
    return llgrid[masked.mask == False]

    def __plot_timings(prefix, show_cb=False, row=0, col=0, the_storage=None):
        _dates = ds["{}_dates.month".format(prefix)][:]

        ax = fig.add_subplot(gs[row, col])

        if the_storage is not None:
            _dates = _dates.where(~np.isnan(the_storage))
            _dates = np.ma.masked_where(the_storage.mask, _dates)

        _dates = maskoceans(lons2d_, lats2d, _dates)

        cs = bmap.pcolormesh(xx, yy, _dates, norm=norm_timings, cmap=cmap_timings)
        cb = bmap.colorbar(cs, location="bottom", format=FuncFormatter(__timing_cb_format_ticklabels))

        if show_cb:
            cb.ax.set_xlabel("month")

            maj_locator = cb.ax.xaxis.get_major_locator()


            print("old tick locs = {}".format(maj_locator.locs))
            maj_locator.locs = __get_new_tick_locs_middle(maj_locator.locs, len(clevs_timings) - 1, shift_direction=-1)
            print("new tick locs = {}".format(maj_locator.locs))


            for tick_line in cb.ax.xaxis.get_ticklines():
                tick_line.set_visible(False)

        cb.ax.set_visible(show_cb)
        ax.set_title("{} timing".format(prefix))
        axes.append(ax)

def plot_only_vegetation_fractions(
        data_path="/RESCUE/skynet3_rech1/huziy/geof_lake_infl_exp/geophys_Quebec_0.1deg_260x260_with_dd_v6_with_ITFS",
        canopy_name="VF", label="QC_10km"):
    r = RPN(data_path)

    veg_fractions = r.get_2D_field_on_all_levels(name=canopy_name)
    print(list(veg_fractions.keys()))

    proj_params = r.get_proj_parameters_for_the_last_read_rec()

    lons, lats = r.get_longitudes_and_latitudes_for_the_last_read_rec()
    print(lons.shape)

    rll = RotatedLatLon(lon1=proj_params["lon1"], lat1=proj_params["lat1"],
                        lon2=proj_params["lon2"], lat2=proj_params["lat2"])

    lon0, lat0 = rll.get_true_pole_coords_in_rotated_system()
    plon, _ = rll.get_north_pole_coords()

    b = Basemap(projection="rotpole", llcrnrlon=lons[0, 0], llcrnrlat=lats[0, 0],
                urcrnrlon=lons[-1, -1], urcrnrlat=lats[-1, -1], lon_0=lon0 - 180,
                o_lon_p=lon0, o_lat_p=lat0)

    lons[lons > 180] -= 360
    for lev in list(veg_fractions.keys()):
        veg_fractions[lev] = maskoceans(lons, lats, veg_fractions[lev], inlands=False)

    x, y = b(lons, lats)
    plot_veg_fractions(x, y, b, veg_fractions, out_image=os.path.join(os.path.dirname(data_path),
                                                                      "veg_fractions_{0}.png".format(label)))

 def __init__(self, filename):
     nc = Dataset(filename)
     self.debug_mode = 1.
     print self.debug_mode
     self.dt = nc.variables["time"][1]-nc.variables["time"][0]
     self.xlon = nc.variables["xlon"][:]
     self.xlat = nc.variables["xlat"][:]
     N_lon = self.xlon.shape[1]
     N_lat = self.xlat.shape[0]
     #self.idx_to_lon = interp1d(numpy.arange(N_lon), self.xlon[0,:])
     #self.idx_to_lat = interp1d(numpy.arange(N_lat), self.xlat[:,0])
     self.idx_to_lon = interp2d(numpy.arange(N_lon), numpy.arange(N_lat), self.xlon)
     self.idx_to_lat = interp2d(numpy.arange(N_lon), numpy.arange(N_lat), self.xlat)
     self.left = self.xlon.min()
     self.right = self.xlon.max()
     self.bottom = self.xlat.min()+18
     self.top = self.xlat.max()-8
     self.m = Basemap(llcrnrlon=self.left,
                      llcrnrlat=self.bottom,
                      urcrnrlon=self.right,
                      urcrnrlat=self.top,
                      projection='cyl', resolution='l')
     self.ocean_mask = maskoceans(self.xlon, self.xlat, self.xlon).mask
     self.time_min = 0.
     self.time_max = 0.
     self.no_ingested = 0
     self.masks = numpy.empty((0,)+self.ocean_mask.shape)
     self.tc_table = numpy.empty((0,5))
     self.time = numpy.empty((0,))

def plot_swe_bfes(runconfig_rea, runconfig_gcm, vname_model="I5", season_to_months=None,
                  bmp_info=None, axes_list=None):

    seasonal_clim_fields_rea = analysis.get_seasonal_climatology_for_runconfig(run_config=runconfig_rea,
                                                                               varname=vname_model, level=0,
                                                                               season_to_months=season_to_months)

    seasonal_clim_fields_gcm = analysis.get_seasonal_climatology_for_runconfig(run_config=runconfig_gcm,
                                                                               varname=vname_model, level=0,
                                                                               season_to_months=season_to_months)

    lons = bmp_info.lons.copy()
    lons[lons > 180] -= 360

    assert len(seasonal_clim_fields_rea) > 0
    season_to_err = OrderedDict()
    for season, field in seasonal_clim_fields_rea.items():
        rea = field
        gcm = seasonal_clim_fields_gcm[season]

        # Mask oceans and lakes

        season_to_err[season] = maskoceans(lons, bmp_info.lats, gcm - rea)
        assert hasattr(season_to_err[season], "mask")

    plot_performance_err_with_cru.plot_seasonal_mean_biases(season_to_error_field=season_to_err,
                                                            varname=vname_model, basemap_info=bmp_info,
                                                            axes_list=axes_list)

def compare_vars(vname_model, vname_to_obs, r_config, season_to_months, bmp_info_agg, axes_list):
    season_to_clim_fields_model = analysis.get_seasonal_climatology_for_runconfig(run_config=r_config,
                                                                                  varname=vname_model, level=0,
                                                                                  season_to_months=season_to_months)

    for season, field in season_to_clim_fields_model.items():
        print(field.shape)
        if vname_model == "PR":
            field *= 1.0e3 * 24 * 3600

    seasonal_clim_fields_obs = vname_to_obs[vname_model]

    lons = bmp_info_agg.lons.copy()
    lons[lons > 180] -= 360

    season_to_err = OrderedDict()
    for season in seasonal_clim_fields_obs:
        season_to_err[season] = season_to_clim_fields_model[season] - seasonal_clim_fields_obs[season]
        season_to_err[season] = maskoceans(lons, bmp_info_agg.lats, season_to_err[season], inlands=False)

    cs = plot_seasonal_mean_biases(season_to_error_field=season_to_err,
                                   varname=vname_model,
                                   basemap_info=bmp_info_agg,
                                   axes_list=axes_list)

    return cs

    def __plot_storage(prefix, show_cb=False, row=0, col=0, plot_deviations=True):

        if plot_deviations:
            clevs = [0, 1e3, 1e4, 1e5, 1.0e6, 1e7, 1e8, 1.0e9]
            clevs = [-c for c in reversed(clevs)][:-1] + clevs
            cmap = cm.get_cmap("bwr", len(clevs) - 1)
        else:
            clevs = [0, 1e3, 1e4, 1e5, 1.0e6, 1e7, 1e8, 1.0e9]
            cmap = cm.get_cmap("YlGnBu", len(clevs) - 1)

        norm = BoundaryNorm(boundaries=clevs, ncolors=len(clevs) - 1)

        _storage = ds["{}_{}".format(prefix, storage_var_name)][:]
        ax = fig.add_subplot(gs[row, col])
        _storage = _storage.where(_storage > storage_lower_limit_m3)
        _storage = maskoceans(lons2d_, lats2d, _storage)
        _storage = np.ma.masked_where(np.isnan(_storage), _storage)


        if plot_deviations:
            cs = bmap.pcolormesh(xx, yy, _storage - bankfull_storage, norm=norm, cmap=cmap)
        else:
            cs = bmap.pcolormesh(xx, yy, _storage, norm=norm, cmap=cmap)


        ext = "both" if plot_deviations else "max"


        cb = bmap.colorbar(cs, location="bottom", format=FuncFormatter(__storage_cb_format_ticklabels), extend=ext)
        cb.ax.set_visible(show_cb)
        cb.ax.set_xlabel(r"${\rm m^3}$")
        ax.set_title(prefix)
        axes.append(ax)
        return _storage

def plotKoppen(Dict):
    """
    Plot a map of the Koppen climate zones used in the correlation analysis.
    """
    m

    empty = np.zeros((27,22))

    empty[23:,14:] = 1.0 #equatorial
    empty[19:23,14:] = 2.0 #tropical
    empty[10:19,16:] = 3.0 #subtropical
    empty[9:17,14:15] = 4.0 #desert

    empty[17:19,14:15] = 5.0#grass1
    empty[17:19,15:16] = 5.0 #grass2
    empty[4:17,15:16] = 5.0 #grass3
    empty[4:9,14:15] = 5.0 #grass4

    empty[4:10,16:] = 6.0 #temperate1
    empty[:4,14:20] = 6.0 #temperate2

    [lonall,latall] = np.meshgrid((Dict['lon']),(Dict['lat']))
    x,y = m(lonall,latall)

    eastern_Aus = maskoceans(lonall,latall,empty)
    eastern_Aus = np.ma.masked_where(eastern_Aus == 0.0,eastern_Aus)
    cs = m.pcolor(x,y,eastern_Aus)
    plt.title("Climate zones in eastern Australia")
    plt.show()
    return eastern_Aus

def plot_field_2d(lons_2d, lats_2d, field_2d, start_lon = -180, end_lon = 0, color_map = None,
                   minmax = (None, None)  ):


    plt.figure()
    m = Basemap(llcrnrlon = start_lon,llcrnrlat = np.min(lats_2d),
                urcrnrlon = end_lon,urcrnrlat = np.max(lats_2d), resolution = 'l')

    m.drawmeridians(range(start_lon,end_lon,10))
    m.drawparallels(range(-90,90,10))


 #   y, x = meshgrid(lats_2d, lons_2d)
#    lons_2d[lons_2d < start_lon] = lons_2d[lons_2d < start_lon] + 360
    x, y = m(lons_2d, lats_2d)


    x -= 360 ###########CONVERTING LONGITUDE TO -180:180
    field_2d = maskoceans(x, y, field_2d)


    m.pcolormesh(x, y, field_2d, cmap = color_map, vmin = minmax[0], vmax = minmax[1])
    m.drawcoastlines()
    #plt.imshow(np.transpose(data[:,:]), origin = 'lower') #for plotting in order to see i,j we supply j,i
    numticks = color_map.N + 1 if color_map != None else 10
    plt.colorbar(ticks = LinearLocator(numticks = numticks), format = '%.01f',
                 orientation = 'vertical', shrink = 0.6)

    def subplot(self, axes, subfigure):
        from mpl_toolkits.basemap import maskoceans
        data, xgrid, ygrid, title, datarange, cbarlabel, map_kwargs = subfigure

        masked_data = maskoceans(xgrid, ygrid, data, inlands=False)
        subfigure = (masked_data, xgrid, ygrid, title, datarange, cbarlab, map_kwargs)
        super(MaskedArrayMapFigure, self).subplot(axes, subfigure)

    def inds_subset(self, lat0=None,lat1=None,lon0=None,lon1=None, maskocean=False, maskland=False):
        ''' return indices of lat,lon arrays within input box '''
        inds=~np.isnan(self.AMF_OMI) # only want non nans
        mlons,mlats=np.meshgrid(self.longitude,self.latitude)
        with np.errstate(invalid='ignore'): # ignore comparisons with NaNs
            if lat0 is not None:
                inds = inds * (mlats >= lat0)
            if lon0 is not None:
                inds = inds * (mlons >= lon0)
            if lat1 is not None:
                inds = inds * (mlats <= lat1)
            if lon1 is not None:
                inds = inds * (mlons <= lon1)

        # true over ocean squares
        oceanmask=maskoceans(mlons,mlats,self.AMF_OMI,inlands=False).mask

        landmask = (~oceanmask)

        # mask ocean if flag is set
        if maskocean:
            inds *= (~oceanmask)
            if __DEBUG__:
                print("oceanmask:")
                print((type(oceanmask),oceanmask.shape))
                print( (inds * (~oceanmask)).shape )
                print((np.sum(oceanmask),np.sum(~oceanmask))) # true for ocean squares!

        if maskland:
            inds *= (~landmask)

        return inds

def draw_map(t_str,
             ax,
             data,
             vmin,
             vmax,
             cmap=plt.get_cmap('Blues'),
             norm=plt.normalize,
             maskoceans_switch=True):

    map = NAMapFigure(t_str=t_str,
                      cb_axis=None,
                      map_axis=ax,
                      fast_or_pretty='pretty',
                      lat_0=49,
                      lon_0=-97,
                      mapwidth=5.8e6,
                      mapheight=5.2e6)

    lon, lat, topo = sp.parse_STEM_coordinates(
        os.path.join(os.getenv('SARIKA_INPUT'), 'TOPO-124x124.nc'))

    if maskoceans_switch:
        data = maskoceans(lon, lat, data, inlands=False, resolution='f')

    cm = map.map.contourf(lon, lat,
                          data,
                          cmap=cmap,
                          latlon=True,
                          norm=norm,
                          vmin=vmin,
                          vmax=vmax)
    return(map, cm)

    def draw_map(self,
                 data,
                 map_axis_idx=None,
                 cb_axis=None,
                 cb_format='%0.2f',
                 cbar_t_str=None,
                 label_lat=False,
                 label_lon=False,
                 vmin=None,
                 vmax=None,
                 midpoint=None,
                 bands_above_mdpt=5,
                 bands_below_mdpt=5,
                 cmap=plt.get_cmap('Blues'),
                 panel_lab='a',
                 extend='neither'):

        """
        returns: stem_pytools.na_map.NAMapFigure object containing the map
        """

        d = STEM_Domain()
        stem_lon = d.get_lon()
        stem_lat = d.get_lat()
        this_cmap, this_norm = midpt_norm.get_discrete_midpt_cmap_norm(
            vmin=vmin,
            vmax=vmax,
            midpoint=midpoint,
            bands_above_mdpt=bands_above_mdpt,
            bands_below_mdpt=bands_below_mdpt,
            extend=extend,
            this_cmap=cmap)

        this_ax = self.ax[map_axis_idx]
        map_obj = na_map.NAMapFigure(map_axis=this_ax,
                                     label_latlon=(label_lat, label_lon),
                                     lon_label_interval=30,
                                     cb_axis=cb_axis,
                                     t_str=None)

        cm = map_obj.map.pcolor(stem_lon, stem_lat,
                                maskoceans(stem_lon, stem_lat, data),
                                cmap=this_cmap,
                                latlon=True,
                                norm=this_norm)
        print 'colorbar format: ' + cb_format
        cbar = this_ax.figure.colorbar(cm,
                                       ax=this_ax,
                                       format=cb_format,
                                       orientation='horizontal')
        if cbar_t_str is not None:
            cbar.ax.set_title(cbar_t_str)
        ticklabs = cbar.ax.get_xticklabels()
        tickvals = np.array([float(x.get_text()) for x in ticklabs])
        cbar.ax.set_xticklabels(tickvals, rotation=-45)
        # place panel label in upper left
        this_ax.text(-0.2, 1.0, panel_lab, transform=this_ax.transAxes)

def compare_swe_diff_for_era40_driven():
    b, lons2d, lats2d = draw_regions.get_basemap_and_coords(llcrnrlat=40.0, llcrnrlon=-145, urcrnrlon=-10)

    lons2d[lons2d > 180] -= 360

    x, y = b(lons2d, lats2d)
    #period
    start_year = 1981
    end_year = 1997
    the_months = [12,1,2]

    levels = [10] + list(range(20, 120, 20)) + [150,200, 300,500,1000]
    cmap = mpl.cm.get_cmap(name="jet_r", lut = len(levels))
    norm = colors.BoundaryNorm(levels, cmap.N)



    swe_obs_manager = SweDataManager(var_name="SWE")
    swe_obs = swe_obs_manager.get_mean(start_year, end_year, months=the_months)
    print("Calculated obs swe")

    swe_obs_interp = swe_obs_manager.interpolate_data_to(swe_obs, lons2d, lats2d)


    axes_list = []

    levels_diff = np.arange(-100, 110, 10)

    #plot model res. (ERA40 driven 1)

    paths = ["data/CORDEX/na/era40_1", "data/CORDEX/na/era40_2"]
    prefixes = ["pmNorthAmerica_0.44deg_ERA40-Int_{0}_1958-1961".format("DJF"),
                "pmNorthAmerica_0.44deg_ERA40-Int2_{0}_1958-1961".format("DJF")
                ]
    pf_kinds = draw_regions.get_permafrost_mask(lons2d, lats2d)
    for i, the_path in enumerate(paths):
        base = os.path.basename(the_path)
        fig = plt.figure()
        ax = plt.gca()
        axes_list.append(ax)


        swe_model_era = CRCMDataManager.get_mean_2d_from_climatologies(path=the_path,
            var_name="I5", file_prefixes=prefixes)
        swe_model_era = maskoceans(lons2d, lats2d, swe_model_era)

        #plot model(ERA40 driven) - obs
        axes_list.append(ax)
        img = b.contourf(x, y, swe_model_era - swe_obs_interp, levels = levels_diff)
        draw_colorbar(fig, img, ax = ax, boundaries=levels_diff)
        ax.set_title("Model ({0} 1958-1961) - Obs.".format(base))

        b.drawcoastlines(ax = ax, linewidth = 0.2)
        b.contour(x, y, pf_kinds, ax = ax, colors = "k")
        fig.savefig("swe_{0}.png".format(base))