Functions

This notebook lists all functions that are defined and used throughout this training material.


Load required libraries

import os
import xarray as xr
from netCDF4 import Dataset
import numpy as np
import glob

# Python libraries for visualization
from matplotlib import pyplot as plt
from matplotlib import pyplot as plt
import matplotlib.colors
from matplotlib.colors import LogNorm
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
import matplotlib.cm as cm

import warnings
warnings.simplefilter(action = "ignore", category = RuntimeWarning)
warnings.simplefilter(action = "ignore", category = FutureWarning)

Data loading and re-shaping functions

generate_xr_from_1d_vec

def generate_xr_from_1d_vec(file, lat_path, lon_path, variable, parameter_name, longname, no_of_dims, unit):
    """ 
    Takes a netCDF4.Dataset or xarray.DataArray object and returns a xarray.DataArray object with latitude / longitude
    information as coordinate information
    
    Parameters:
        file(netCDF4 data file or xarray.Dataset): AC SAF or IASI Level 2 data file, loaded a netCDF4.Dataset or xarray.DataArray
        lat_path(str): internal path of the data file to the latitude information, e.g. 'GEOLOCATION/LatitudeCentre'
        lon_path(str): internal path of the data file to the longitude information, e.g. 'GEOLOCATION/LongitudeCentre'
        variable(array): extracted variable of interested
        parameter_name(str): parameter name, preferably extracted from the data file
        longname(str): Long name of the parameter, preferably extracted from the data file
        no_of_dims(int): Define the number of dimensions of your input array
        unit(str): Unit of the parameter, preferably extracted from the data file
    
    Returns:
        1 or 2-dimensional (depending on the given number of dimensions) xarray.DataArray  with latitude / longitude information 
        as coordinate information
    """
    latitude = file[lat_path]
    longitude = file[lon_path]
    param = variable 
    
    if (no_of_dims==1):
        param_da = xr.DataArray(
            param[:],
            dims=('ground_pixel'),
            coords={
                'latitude': ('ground_pixel', latitude.data),
                'longitude': ('ground_pixel', longitude.data)
            },
            attrs={'long_name': longname, 'units': unit},
            name=parameter_name
        )
    else:
        param_da = xr.DataArray(
            param[:],
            dims=["x","y"],
            coords={
                'latitude':(['x','y'],latitude[:]),
                'longitude':(['x','y'],longitude[:])
            },
            attrs={'long_name': longname, 'units': unit},
            name=parameter_name
        )
        
    return param_da

load_l2_data_xr

def load_l2_data_xr(directory, internal_filepath, parameter, lat_path, lon_path, no_of_dims, 
                    paramname, unit, longname):
    """ 
    Loads a Metop-A/B Level 2 dataset in HDF format and returns a xarray.DataArray with all the ground pixels of all directory 
    files. Uses function 'generate_xr_from_1d_vec' to generate the xarray.DataArray.
    
    Parameters:
        directory(str): directory where the HDF files are stored
        internal_filepath(str): internal path of the data file that is of interest, e.g. TOTAL_COLUMNS
        parameter(str): paramter that is of interest, e.g. NO2
        lat_path(str): name of latitude variable
        lon_path(str): name of longitude variable
        no_of_dims(int): number of dimensions of input array
        paramname(str): name of parameter
        unit(str): unit of the parameter, preferably taken from the data file
        longname(str): longname of the parameter, preferably taken from the data file
    
    Returns:
        1 or 2-dimensional xarray.DataArray with latitude / longitude information as coordinate information
    """
    filelist = glob.glob(directory+'/*')
    datasets = []

    for i in filelist:
        tmp=Dataset(i)
        param=tmp[internal_filepath+'/'+parameter]
        da_tmp= generate_xr_from_1d_vec(tmp,lat_path, lon_path,
                                param, paramname, longname, no_of_dims, unit)
        if(no_of_dims==1):
            datasets.append(da_tmp)
        else:
            da_tmp_st = da_tmp.stack(ground_pixel=('x','y'))
            datasets.append(da_tmp_st)

    return xr.concat(datasets, dim='ground_pixel')

generate_geographical_subset

def generate_geographical_subset(xarray, latmin, latmax, lonmin, lonmax, reassign=False):
    """ 
    Generates a geographical subset of a xarray.DataArray and if kwarg reassign=True, shifts the longitude grid 
    from a 0-360 to a -180 to 180 deg grid.
    
    Parameters:
        xarray(xarray.DataArray): a xarray DataArray with latitude and longitude coordinates
        latmin, latmax, lonmin, lonmax(int): lat/lon boundaries of the geographical subset
        reassign(boolean): default is False
        
    Returns:
        Geographical subset of a xarray.DataArray.
    """   
    if(reassign):
        xarray = xarray.assign_coords(longitude=(((xarray.longitude + 180) % 360) - 180))
    return xarray.where((xarray.latitude < latmax) & (xarray.latitude > latmin) & (xarray.longitude < lonmax) & (xarray.longitude > lonmin),drop=True)

generate_masked_array

def generate_masked_array(xarray, mask, threshold, operator, drop=True):
    """ 
    Applies a mask (e.g. a cloud mask) onto a given xarray.DataArray, based on a given threshold and operator.
    
    Parameters:
        xarray(xarray DataArray): a three-dimensional xarray.DataArray object
        mask(xarray DataArray): 1-dimensional xarray.DataArray, e.g. cloud fraction values
        threshold(float): any number specifying the threshold
        operator(str): operator how to mask the array, e.g. '<', '>' or '!='
        drop(boolean): default is True
        
    Returns:
        Masked xarray.DataArray with NaN values dropped, if kwarg drop equals True
    """
    if(operator=='<'):
        cloud_mask = xr.where(mask < threshold, 1, 0) #Generate cloud mask with value 1 for the pixels we want to keep
    elif(operator=='!='):
        cloud_mask = xr.where(mask != threshold, 1, 0)
    elif(operator=='>'):
        cloud_mask = xr.where(mask > threshold, 1, 0)
    else:
        cloud_mask = xr.where(mask == threshold, 1, 0)
            
    xarray_masked = xr.where(cloud_mask ==1, xarray, np.nan) #Apply mask onto the DataArray
    xarray_masked.attrs = xarray.attrs #Set DataArray attributes 
    if(drop):
        return xarray_masked[~np.isnan(xarray_masked)] #Return masked DataArray
    else:
        return xarray_masked

load_masked_l2_da

def load_masked_l2_da(directory, internal_filepath, parameter, lat_path, lon_path, no_of_dims, 
                      paramname, longname, unit, threshold, operator):
    """ 
    Loads a Metop-A/B Gome-2 Level 2 data and cloud fraction information and 
    returns a masked xarray.DataArray. It combines the functions `load_l2_data_xr` and `generate_masked_array`.
    
    Parameters:
        directory(str): Path to directory with Level 2 data files.
        internal_filepath(str): Internal file path under which the parameters are strored, e.g. TOTAL_COLUMNS
        parameter(str): atmospheric parameter, e.g. NO2
        lat_path(str): name of the latitude variable within the file
        lon_path(str): path to the longitude variable within the file
        no_of_dims(int): specify the number of dimensions, 1 or 2
        paramname(str): parameter name
        longname(str): long name of the parameter that shall be used
        unit(str): unit of the parameter
        threshold(float): any number specifying the threshold
        operator(str): operator how to mask the xarray.DataArray, e.g. '<', '>' or '!='
        
    Returns:
        Masked xarray.DataArray keeping NaN values (drop=False)
    """  
    da = load_l2_data_xr(directory, 
                         internal_filepath, 
                         parameter, 
                         lat_path, 
                         lon_path, 
                         no_of_dims, 
                         paramname, 
                         unit, 
                         longname)
    
    cloud_fraction = load_l2_data_xr(directory, 
                                     'CLOUD_PROPERTIES', 
                                     'CloudFraction', 
                                     lat_path, 
                                     lon_path, 
                                     no_of_dims, 
                                     'CloudFraction', 
                                     unit='-', 
                                     longname='Cloud Fraction')
    
    return generate_masked_array(da, cloud_fraction, threshold, operator, drop=False)

select_channels_for_rgb

def select_channels_for_rgb(xarray, red_channel, green_channel, blue_channel):
    """ 
    Selects the channels / bands of a multi-dimensional xarray for red, green and blue composite based on Sentinel-3
    OLCI Level 1B data.
    
    Parameters:
        xarray(xarray.Dataset): xarray.Dataset object that stores the different channels / bands.
        red_channel(str): Name of red channel to be selected
        green_channel(str): Name of green channel to be selected
        blue_channel(str): Name of blue channel to be selected

    Returns:
        Three xarray DataArray objects with selected channels / bands
    """  
    return xarray[red_channel], xarray[green_channel], xarray[blue_channel]

normalize

def normalize(array):
    """ 
    Normalizes a numpy array / xarray.DataArray object to values between 0 and 1.
    
    Parameters:
        xarray(numpy array or xarray.DataArray): xarray.DataArray or numpy array object whose values should be
        normalized.

    Returns:
        xarray.DataArray with normalized values
    """ 
    array_min, array_max = array.min(), array.max()
    return ((array - array_min)/(array_max - array_min))

slstr_frp_gridding

def slstr_frp_gridding(parameter_array, parameter, lat_min, lat_max, lon_min, lon_max, 
                       sampling_lat_frp_grid, sampling_lon_frp_grid, n_fire, lat_frp, lon_frp, **kwargs):
    """ 
    Produces gridded data of Sentinel-3 SLSTR NRT Fire Radiative Power Data
    
    Parameters:
        parameter_array(xarray.DataArray): xarray.DataArray with extracted data variable of fire occurences
        parameter(str): NRT S3 frp channel - either `mwir`, `swir` or `swir_nosaa`
        lat_min, lat_max, lon_min, lon_max(float): Floats of geographical bounding box
        sampling_lat_frp_grid, sampling_long_frp_grid(float): Float of grid cell size
        n_fire(int): Number of fire occurences
        lat_frp(xarray.DataArray): Latitude values of occurred fire events
        lon_frp(xarray.DataArray): Longitude values of occurred fire events
        **kwargs: additional keyword arguments to be added. Required for parameter `swir_nosaa`, where the function
                  requires the xarray.DataArray with the SAA FLAG information.  

    Returns:
        the gridded xarray.Data Array and latitude and longitude grid information
    """ 
    n_lat = int( (np.float32(lat_max) - np.float32(lat_min)) / sampling_lat_frp_grid ) + 1 # Number of rows per latitude sampling
    n_lon = int( (np.float32(lon_max) - np.float32(lon_min)) / sampling_lon_frp_grid ) + 1 # Number of lines per longitude sampling

    
    slstr_frp_gridded = np.zeros( [n_lat, n_lon], dtype='float32' ) - 9999.

    lat_grid = np.zeros( [n_lat, n_lon], dtype='float32' ) - 9999.
    lon_grid = np.zeros( [n_lat, n_lon], dtype='float32' ) - 9999.
    
    if (n_fire >= 0):
    
    # Loop on i_lat: begins
        for i_lat in range(n_lat):
                    
        # Loop on i_lon: begins
            for i_lon in range(n_lon):
                        
                lat_grid[i_lat, i_lon] = lat_min + np.float32(i_lat) * sampling_lat_frp_grid + sampling_lat_frp_grid / 2.
                lon_grid[i_lat, i_lon] = lon_min + np.float32(i_lon) * sampling_lon_frp_grid + sampling_lon_frp_grid / 2.
                            
            # Gridded SLSTR frp MWIR Night - All days
                if(parameter=='swir_nosaa'):
                    flag_frp_swir_saa_nc = kwargs.get('flag', None)
                    mask_grid = np.where( 
                        (lat_frp[:] >= lat_min + np.float32(i_lat) * sampling_lat_frp_grid)  & 
                        (lat_frp[:] < lat_min + np.float32(i_lat+1) * sampling_lat_frp_grid) & 
                        (lon_frp[:] >= lon_min + np.float32(i_lon) * sampling_lon_frp_grid)  & 
                        (lon_frp[:] < lon_min + np.float32(i_lon+1) * sampling_lon_frp_grid) &
                        (parameter_array[:] != -1.) & (flag_frp_swir_saa_nc[:] == 0), False, True)
                else:
                    mask_grid = np.where( 
                        (lat_frp[:] >= lat_min + np.float32(i_lat) * sampling_lat_frp_grid)  & 
                        (lat_frp[:] < lat_min + np.float32(i_lat+1) * sampling_lat_frp_grid) & 
                        (lon_frp[:] >= lon_min + np.float32(i_lon) * sampling_lon_frp_grid)  & 
                        (lon_frp[:] < lon_min + np.float32(i_lon+1) * sampling_lon_frp_grid) &
                        (parameter_array[:] != -1.),  False, True)
                            
                masked_slstr_frp_grid = np.ma.array(parameter_array[:], mask=mask_grid)
                            
                if len(masked_slstr_frp_grid.compressed()) != 0:
                    slstr_frp_gridded[i_lat, i_lon]  = np.sum(masked_slstr_frp_grid.compressed())
    return slstr_frp_gridded, lat_grid, lon_grid
    

df_subset

def df_subset(df,low_bound1, high_bound1, low_bound2, high_bound2):
    """ 
    Creates a subset of a pandas.DataFrame object with time-series information
    
    Parameters:
        df(pandas.DataFrame): pandas.DataFrame with time-series information
        low_bound1(str): dateTime string, e.g. '2018-11-30'
        high_bound1(str): dateTime string, e.g. '2018-12-01'
        low_bound2(str): dateTime string, e.g. '2019-12-30'
        high_bound2(str): dateTime string, e.g. '2020-01-15'

    Returns:
        the subsetted time-series as pandas.DataFrame object
    """ 
    return df[(df.index>low_bound1) & (df.index<high_bound1)], df[(df.index>low_bound2) & (df.index<high_bound2)]

Data visualization functions

visualize_imshow

def visualize_imshow(data_array, projection, conversion_factor, color_scale, vmin, vmax, lonmin, lonmax, latmin, latmax, unit, set_global=False, log_scale=False):
    """ 
    Visualizes a numpy MaskedArray with matplotlib's 'imshow' function.
    
    Parameters:
        data_array (numpy MaskedArray): any numpy MaskedArray, e.g. loaded with the NetCDF library and the Dataset function
        projection (str): a projection provided by the cartopy library, e.g. ccrs.PlateCarree()
        conversion_factor (float): any number to convert the DataArray values
        color_scale(str): string taken from matplotlib's color ramp reference  
        vmin (int): minimum number on visualisation legend
        vmax (int): maximum number on visualisation legend
        lonmin, lonmax, latmin, latmax (float): geographic boundary values 
        unit (str): define unit of the plot to be added to the colorbar
        set_global (logical): set True, if the plot shall have a global coverage
        log_scale (logical): set True, if the color_scale shall have a logarithmic scaling
    """
    fig=plt.figure(figsize=(20, 12))

    ax=plt.axes(projection=projection)
    ax.coastlines()

    if(set_global):
        ax.set_global()
        ax.gridlines()
    
    if (projection==ccrs.PlateCarree()):
        ax.set_extent([lonmin, lonmax, latmin, latmax], projection)
        gl = ax.gridlines(draw_labels=True, linestyle='--')
        gl.xlabels_top=False
        gl.ylabels_right=False
        gl.xformatter=LONGITUDE_FORMATTER
        gl.yformatter=LATITUDE_FORMATTER
        gl.xlabel_style={'size':14}
        gl.ylabel_style={'size':14}

    if(log_scale):
        img1 = plt.imshow(data_array[:]*conversion_factor,
                          cmap=color_scale,
                          aspect='auto',
                          norm=matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax))
    else:
        img1 = plt.imshow(data_array[:]*conversion_factor,
                          cmap=color_scale,
                          vmin=vmin,
                          vmax=vmax,
                          aspect='auto')

    cbar = fig.colorbar(img1, ax=ax, orientation='horizontal', fraction=0.04, pad=0.1)
    cbar.set_label(str(conversion_factor) + ' ' + unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)
    
    plt.show()

visualize_scatter

def visualize_scatter(xr_dataarray, conversion_factor, projection, vmin, vmax, point_size, color_scale, unit, 
                      title):
    """ 
    Visualizes a xarray.DataArray in a given projection using matplotlib's scatter function.
    
    Parameters:
        xr_dataarray(xarray.DataArray): a one-dimensional xarray DataArray object with latitude and longitude information as coordinates
        conversion_factor(int): any number to convert the DataArray values
        projection(str): choose one of cartopy's projection, e.g. ccrs.PlateCarree()
        vmin(int): minimum number on visualisation legend
        vmax(int): maximum number on visualisation legend
        point_size(int): size of marker, e.g. 5
        color_scale(str): string taken from matplotlib's color ramp reference
        unit(str): define the unit to be added to the color bar
        title(str): define title of the plot
    """
    fig = plt.figure(figsize=(16, 10))
    ax = fig.add_subplot(1, 1, 1, projection=projection)

    ax.coastlines()
    
    if (projection==ccrs.PlateCarree()):
        gl = ax.gridlines(draw_labels=True, linestyle='--')
        gl.top_labels=False
        gl.right_labels=False
        gl.xformatter=LONGITUDE_FORMATTER
        gl.yformatter=LATITUDE_FORMATTER
        gl.xlabel_style={'size':14}
        gl.ylabel_style={'size':14}

    # plot pixel positions
    img = ax.scatter(
        xr_dataarray.longitude.data,
        xr_dataarray.latitude.data,
        c=xr_dataarray.data*conversion_factor,
        cmap=plt.cm.get_cmap(color_scale),
        marker='o',
        s=point_size,
        transform=ccrs.PlateCarree(),
        vmin=vmin,
        vmax=vmax
    )

    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    plt.xlabel("Longitude", fontsize=16)
    plt.ylabel("Latitude", fontsize=16)
    cbar = fig.colorbar(img, ax=ax, orientation='horizontal', fraction=0.04, pad=0.1)
    cbar.set_label(unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)
    ax.set_title(title, fontsize=20, pad=20.0)
    plt.show()

visualize_pcolormesh

def visualize_pcolormesh(data_array, longitude, latitude, projection, color_scale, unit, long_name, vmin, vmax, 
                         set_global=True, lonmin=-180, lonmax=180, latmin=-90, latmax=90):
    """ 
    Visualizes a xarray.DataArray with matplotlib's pcolormesh function.
    
    Parameters:
        data_array(xarray.DataArray): xarray.DataArray holding the data values
        longitude(xarray.DataArray): xarray.DataArray holding the longitude values
        latitude(xarray.DataArray): xarray.DataArray holding the latitude values
        projection(str): a projection provided by the cartopy library, e.g. ccrs.PlateCarree()
        color_scale(str): string taken from matplotlib's color ramp reference
        unit(str): the unit of the parameter, taken from the NetCDF file if possible
        long_name(str): long name of the parameter, taken from the NetCDF file if possible
        vmin(int): minimum number on visualisation legend
        vmax(int): maximum number on visualisation legend
        set_global(boolean): optional kwarg, default is True
        lonmin,lonmax,latmin,latmax(float): optional kwarg, set geographic extent is set_global kwarg is set to 
                                            False

    """
    fig=plt.figure(figsize=(20, 10))

    ax = plt.axes(projection=projection)
   
    img = ax.pcolormesh(longitude, latitude, data_array, 
                        cmap=plt.get_cmap(color_scale), transform=ccrs.PlateCarree(),
                        vmin=vmin,
                        vmax=vmax,
                        shading='auto')

    ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=1)
    ax.add_feature(cfeature.COASTLINE, edgecolor='black', linewidth=1)

    if (projection==ccrs.PlateCarree()):
        ax.set_extent([lonmin, lonmax, latmin, latmax], projection)
        gl = ax.gridlines(draw_labels=True, linestyle='--')
        gl.top_labels=False
        gl.right_labels=False
        gl.xformatter=LONGITUDE_FORMATTER
        gl.yformatter=LATITUDE_FORMATTER
        gl.xlabel_style={'size':14}
        gl.ylabel_style={'size':14}

    if(set_global):
        ax.set_global()
        ax.gridlines()

    cbar = fig.colorbar(img, ax=ax, orientation='horizontal', fraction=0.04, pad=0.1)
    cbar.set_label(unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)
    ax.set_title(long_name, fontsize=20, pad=20.0)

 #   plt.show()
    return fig, ax

visualize_s3_pcolormesh

def visualize_s3_pcolormesh(color_array, array, latitude, longitude, title):
    """ 
    Visualizes a xarray.DataArray or numpy.MaskedArray (Sentinel-3 OLCI Level 1 data) with matplotlib's pcolormesh function as RGB image.
    
    Parameters:
        color_array (numpy.MaskedArray): any numpy.MaskedArray, e.g. loaded with the NetCDF library and the Dataset function
        array(numpy.Array): numpy.Array to get dimensions of the resulting plot
        longitude (numpy.Array): array with longitude values
        latitude (numpy.Array) : array with latitude values
        title (str): title of the resulting plot
    """
    fig=plt.figure(figsize=(20, 12))

    ax=plt.axes(projection=ccrs.Mercator())
    ax.coastlines()

    gl = ax.gridlines(draw_labels=True, linestyle='--')
    gl.top_labels=False
    gl.right_labels=False
    gl.xformatter=LONGITUDE_FORMATTER
    gl.yformatter=LATITUDE_FORMATTER
    gl.xlabel_style={'size':14}
    gl.ylabel_style={'size':14}

    img1 = plt.pcolormesh(longitude, latitude, array*np.nan, color=color_array,
                          clip_on = True,
                          edgecolors=None,
                          zorder=0,
                          transform=ccrs.PlateCarree())
    ax.set_title(title, fontsize=20, pad=20.0)
    plt.show()

visualize_s3_frp

def visualize_s3_frp(data, lat, lon, unit, longname, textstr_1, textstr_2, vmax):
    """ 
    Visualizes a numpy.Array (Sentinel-3 SLSTR NRT frp data) with matplotlib's pcolormesh function and adds two
    text boxes to the plot.
    
    Parameters:
        data(numpy.MaskedArray): any numpy MaskedArray, e.g. loaded with the NetCDF library and the Dataset function
        lat(numpy.Array): array with longitude values
        lon(numpy.Array) : array with latitude values
        unit(str): unit of the resulting plot
        longname(str): Longname to be used as title
        textstr_1(str): String to fill box 1
        textstr_2(str): String to fill box 2
        vmax(float): Maximum value of color scale
    """
    fig=plt.figure(figsize=(20, 15))

    ax = plt.axes(projection=ccrs.PlateCarree())

    img = plt.pcolormesh(lon, lat, data, 
                        cmap=cm.autumn_r, transform=ccrs.PlateCarree(),
                        vmin=0,
                        vmax=vmax)

    ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=1)
    ax.add_feature(cfeature.COASTLINE, edgecolor='black', linewidth=1)

    gl = ax.gridlines(draw_labels=True, linestyle='--')
    gl.bottom_labels=False
    gl.right_labels=False
    gl.xformatter=LONGITUDE_FORMATTER
    gl.yformatter=LATITUDE_FORMATTER
    gl.xlabel_style={'size':14}
    gl.ylabel_style={'size':14}

    cbar = fig.colorbar(img, ax=ax, orientation='horizontal', fraction=0.029, pad=0.025)
    cbar.set_label(unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)
    ax.set_title(longname, fontsize=20, pad=40.0) 

    props = dict(boxstyle='square', facecolor='white', alpha=0.5)

    # place a text box on the right side of the plot
    ax.text(1.1, 0.9, textstr_1, transform=ax.transAxes, fontsize=16,
        verticalalignment='top', bbox=props)

    props = dict(boxstyle='square', facecolor='white', alpha=0.5)

    # place a text box in upper left in axes coords
    ax.text(1.1, 0.85, textstr_2, transform=ax.transAxes, fontsize=16,
            verticalalignment='top', bbox=props)
    plt.show()

visualize_s3_aod

def visualize_s3_aod(aod_ocean, aod_land, latitude, longitude, title, unit, vmin, vmax, color_scale, projection):
    """ 
    Visualizes two xarray.DataArrays from the Sentinel-3 SLSTR NRT AOD dataset onto the same plot with 
    matplotlib's pcolormesh function.
    
    Parameters:
        aod_ocean(xarray.DataArray): xarray.DataArray with the Aerosol Optical Depth for ocean values
        aod_land(xarray.DataArray): xarray.DataArray with Aerosol Optical Depth for land values
        longitude(xarray.DataArray): xarray.DataArray holding the longitude values
        latitude(xarray.DataArray): xarray.DataArray holding the latitude values
        title(str): title of the resulting plot
        unit(str): unit of the resulting plot
        vmin(int): minimum number on visualisation legend
        vmax(int): maximum number on visualisation legend
        color_scale(str): string taken from matplotlib's color ramp reference
        projection(str): a projection provided by the cartopy library, e.g. ccrs.PlateCarree()
    """
    fig=plt.figure(figsize=(12, 12))

    ax=plt.axes(projection=projection)
    ax.coastlines(linewidth=1.5, linestyle='solid', color='k', zorder=10)

    gl = ax.gridlines(draw_labels=True, linestyle='--')
    gl.top_labels=False
    gl.right_labels=False
    gl.xformatter=LONGITUDE_FORMATTER
    gl.yformatter=LATITUDE_FORMATTER
    gl.xlabel_style={'size':12}
    gl.ylabel_style={'size':12}


    img1 = plt.pcolormesh(longitude, latitude, aod_ocean, transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax, cmap=color_scale)
    img2 = plt.pcolormesh(longitude, latitude, aod_land, transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax, cmap=color_scale)
    ax.set_title(title, fontsize=20, pad=20.0)

    cbar = fig.colorbar(img1, ax=ax, orientation='vertical', fraction=0.04, pad=0.05)
    cbar.set_label(unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)

    plt.show()

References