<jupyter_start><jupyter_text>OLCI spatial plotting, quality control and data interrogation     Version: 3.0    Date:    10/07/2020    Author:  <NAME> (InnoFlair / Plymouth Marine Laboratory) and <NAME> (EUMETSAT)    Credit:  This code was developed for EUMETSAT under contracts for the European Commission Copernicus              programme.    License: This code is offered as open source and free-to-use in the public domain,              with no warranty, under the MIT license associated with this code repository. This aim of this code is to introduce you to Python and to simply import a netCDF file in to your python workspace, conduct some basic operations, and plot an image. In this case, we will be using a level-2 OLCI image, but the script can be easily adapted to plot any netCDF variable.    Get the WEkEO User credentialsIf you want to download the data to use this notebook, you will need WEkEO User credentials. If you do not have these, you can register here. The first step in any python code is usually to import libraries that you will need. Libraries are usually code modules that perform specific tasks or provide specific capability (e.g. statistical analysis or plotting routines). In this case we will import the xarray library for handling netCDF files, the numpy library which will help to conduct various operations on the data, and the matplotlib plotting library to generate some images. We will also import the os library, that allows python access to some command-line-eqsue capability like 'list directory', as well as the python library that governs the reporting of warning (so that we can turn them off here, and make the code run without being so 'noisy').<jupyter_code>%matplotlib inline

# libraries are imported here, and we can import any library with an alias that allows us easy access to them later.
import sys
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import os
import glob
from IPython.core.display import display, HTML
import json
from zipfile import ZipFile

import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>We are also going to import the WEkEO Harmonised Data Access API, which will allow us to get remote data if required.<jupyter_code>sys.path.append(os.path.join(os.path.dirname(os.path.dirname(os.getcwd())),'wekeo-hda'))
import hda_api_functions as hapi<jupyter_output><empty_output><jupyter_text>Usually we also define functions at the top of a Python script. Functions are routines that can be called elsewhere in our script and perform a specific task. Typically we would use a function to take care of any process that we are going to perform more than once. The box below defines a function that will mask our data according to quality flags. We will call this function later on.<jupyter_code>def flag_data_fast(flags_we_want, flag_names, flag_values, flag_data, flag_type='WQSF'):
    flag_bits = np.uint64()
    if flag_type == 'SST':
        flag_bits = np.uint8()
    elif flag_type == 'WQSF_lsb':
        flag_bits = np.uint32()
    
    for flag in flags_we_want:
        try:
            flag_bits = flag_bits | flag_values[flag_names.index(flag)]
        except:
            print(flag + " not present")
    
    return (flag_data & flag_bits) > 0<jupyter_output><empty_output><jupyter_text>Now we will start our script, proper.To run this script, we will download some data from WEkEO harmonised data access. WEkEO provides access to a huge number of datasets through its **'harmonised-data-access'** API. This allows us to query the full data catalogue and download data quickly and directly onto the Jupyter Lab. You can search for what data is available hereIn order to use the HDA-API we need to provide some authentication credentials, which comes in the form of an API key and API token. In this notebook we have provided functions so you can retrieve the API key and token you need directly. You can find out more about this process in the notebook on HDA access (wekeo_harmonized_data_access_api.ipynb) that can be found in the **wekeo-hda** folder on your Jupyterlab.We will also define a few other parameters including where to download the data to, and if we want the HDA-API functions to be verbose. **Lastly, we will tell the notebook where to find the query we will use to find the data.** These 'JSON' queries are what we use to ask WEkEO for data. They have a very specific form, but allow us quite fine grained control over what data to get. You can find the example one that we will use here: **JSON_templates/EO_EUM_DAT_SENTINEL-3_OL_2_WRR___.json**<jupyter_code># your WEkEO API username and password (needs to be in '  ')
user_name = 'USERNAME'
password = 'PASSWORD'

# Generate an API key
api_key = hapi.generate_api_key(user_name, password)
display(HTML('Your API key is: <b>'+api_key+'</b>'))

# set this key to true to download data.
download_data = True
# where the data should be downloaded to:
download_dir_path = os.path.join(os.getcwd(),'products')
# where we can find our data query form:
JSON_query_dir = os.path.join(os.getcwd(),'JSON_templates')
# HDA-API loud and noisy?
verbose = False

# make the output directory if required
if not os.path.exists(download_dir_path):
    os.makedirs(download_dir_path)<jupyter_output><empty_output><jupyter_text>Now we have set how we want the script to run, we are ready to get some data. We start this process by telling the script what kind of data we want. In this case, this is OLCI L2 reduced resolution data, which has the following designation on WEkEO: **EO:EUM:DAT:SENTINEL-3:OL_2_WRR___**.<jupyter_code># OLCI REDUCED RESOLUTION LEVEL 2 DATASET ID
dataset_id = "EO:EUM:DAT:SENTINEL-3:OL_2_WRR___"<jupyter_output><empty_output><jupyter_text>Here, we use this dataset_id to find the correct, locally stored JSON query file which describes the data we want. The query file is called: **JSON_templates/cyano/EO_EUM_DAT_SENTINEL-3_OL_2_WRR___.json**You can edit this query if you want to get different data, but be aware of asking for too much data - you could be here a while and might run out of space to use this data in the JupyterLab. The box below gets the correct query file.<jupyter_code># find query file
JSON_query_file = os.path.join(JSON_query_dir,dataset_id.replace(':','_')+".json")
if not os.path.exists(JSON_query_file):
    print('Query file ' + JSON_query_file + ' does not exist')
else:
    print('Found JSON query file for '+dataset_id)<jupyter_output><empty_output><jupyter_text>Now we have a query, we need to launch it to WEkEO to get our data. The box below takes care of this through the following steps:    1. initialise our HDA-API    2. get an access token for our data    3. accepts the WEkEO terms and conditions    4. loads our JSON query into memory    5. launches our search    6. waits for our search results    7. gets our result list    8. downloads our dataThis is quite a complex process, so much of the functionality has been buried 'behind the scenes'. If you want more information, you can check out the **wekeo-hda** tool kit in the parent training directory. The code below will report some information as it runs. At the end, it should tell you that one product has been downloaded.<jupyter_code>if download_data:
    HAPI_dict = hapi.init(dataset_id, api_key, download_dir_path)
    HAPI_dict = hapi.get_access_token(HAPI_dict)
    HAPI_dict = hapi.acceptTandC(HAPI_dict)

    # load the query
    with open(JSON_query_file, 'r') as f:
        query = json.load(f)

    # launch job
    print('Launching job...')
    HAPI_dict = hapi.get_job_id(HAPI_dict, query)

    # check results
    print('Getting results...')
    HAPI_dict = hapi.get_results_list(HAPI_dict)
    HAPI_dict = hapi.get_order_ids(HAPI_dict)

    # download data
    print('Downloading data...')
    HAPI_dict = hapi.download_data(HAPI_dict, file_extension='.zip')<jupyter_output><empty_output><jupyter_text>Sentinel data is usually distributed as a zip file, which contains the SAFE format data within. To use this, we must unzip the file. The box below handles this.<jupyter_code>if download_data:
    # unzip file
    for filename in HAPI_dict['filenames']:
        if os.path.splitext(filename)[-1] == '.zip':
            print('Unzipping file')
            try:
                with ZipFile(filename, 'r') as zipObj:
                    # Extract all the contents of zip file in current directory
                    zipObj.extractall(os.path.dirname(filename))

                # clear up the zip file
                os.remove(filename)
            except:
                print("Failed to unzip....")
if download_data:
    unzipped_file = HAPI_dict['filenames'][0].replace('.zip','.SEN3')
else:
    unzipped_file = glob.glob(os.path.join(download_dir_path,'*OL_2_WRR*.SEN3'))

input_root    = os.path.dirname(unzipped_file)
input_path    = os.path.basename(unzipped_file)
file_name_chl = 'chl_nn.nc'<jupyter_output><empty_output><jupyter_text>We'll quickly check, in the next box, if your data path is ok, and that the data file exists check.<jupyter_code># quick path length check (some windows versions have a problem with long file paths)
if len(os.path.join(input_root,input_path,file_name_chl)) > 259 \
  or len(os.path.join(input_root,input_path,file_name_chl)) > 248:
    print('Beware, your path name is quite long. Consider moving your data to a new directory')
else:
    print('Path length name seems fine')
    
if os.path.exists(os.path.join(input_root,input_path,file_name_chl)):
    print('Found the required data file')
else:
    print('Data file missing. Please check your path and file name')<jupyter_output><empty_output><jupyter_text>We read the file using functions from the netCDF4 (alias "nc") library. Note that to use a library in python you use the imported alias followed by a dot, and then the function you want (e.g. nc.Dataset).<jupyter_code>OLCI_file = xr.open_dataset(os.path.join(input_root,input_path,file_name_chl))<jupyter_output><empty_output><jupyter_text>To access a variable you can use the following command, where the name of the variable you are interested in, follows the hash. If you remove the hash in the following box, put the cursor after the dot and hit 'tab' you will be presented with a list of all of the variables and methods associated with the OLCI_file object. Python is an 'object orientatated' language, which means that all objects have relevant methods associated with them. note: If you want to run all this code in one go, remember to put the hash back at the start of this line beforehand.<jupyter_code>#OLCI_file.<jupyter_output><empty_output><jupyter_text>So, lets load in some data, and then close our data file<jupyter_code>CHL     = OLCI_file.CHL_NN.data
OLCI_file.close()<jupyter_output><empty_output><jupyter_text>You can look at the variables in your workspace in interactive python environments (like this, or ipython) by typing 'whos'. This will tell you the name of the variable, it's type and then information on it, such as its size and shape.<jupyter_code>#whos<jupyter_output><empty_output><jupyter_text>Lets take a look at our data..<jupyter_code>plt.imshow(CHL);<jupyter_output><empty_output><jupyter_text>This is not the prettiest plot - Python can do much better. For a start, we may wish to look at a smaller area. We'll do this now, using the relevant indexes for area of data you wish to use.<jupyter_code>row1=3000
row2=7000
col1=0
col2=3000
grid_factor=2
CHL_subset = CHL[row1:row2:grid_factor, col1:col2:grid_factor]
plt.imshow(CHL_subset);<jupyter_output><empty_output><jupyter_text>You will notice a few problems with displaying plots like this. Firstly - they don't look very pretty (the colour scheme is not ideal, it is hard to see the coastline, and you can't differentiate the land from cloud), and secondly - the axes don't provide any information on the location (other than within the array) and there is no colour bar. To make a better plot we will need to add a few more tools to the libraries/modules we've currently imported. Below are a few lines of code to import cartopy (which will help us make a better, map based plot) and a few other tools to tweak how the data is displayed. The Cartopy is module allows us to use map projections to display data in a geographically relevant way. For those that are familiar to python, Cartopy has largely replaced the Basemap library.<jupyter_code>import cartopy.crs as ccrs
import cartopy.feature as cfeature

land_resolution = '50m'
land_poly = cfeature.NaturalEarthFeature('physical', 'land', land_resolution,
                                        edgecolor='k',
                                        facecolor=cfeature.COLORS['land'])<jupyter_output><empty_output><jupyter_text>We will also need to load other data to make the plot - the longitude and latitude data associated with each pixel of the chlorophyll data. This data can be found in the geo_coordinates.nc file, within each S3 OLCI L2 folder. We load this in a very similar way to how we loaded the chlorophyll data, just with different file and variable names. The data path remains the same, referring the folder that contains all the netcdf files.<jupyter_code>file_name_geo = 'geo_coordinates.nc'
GEO_file      = xr.open_dataset(os.path.join(input_root,input_path,file_name_geo))
LAT           = GEO_file.variables['latitude'][:]
LON           = GEO_file.variables['longitude'][:]
LAT_subset    = LAT[row1:row2:grid_factor, col1:col2:grid_factor]
LON_subset    = LON[row1:row2:grid_factor, col1:col2:grid_factor]
GEO_file.close()<jupyter_output><empty_output><jupyter_text>Then we need to initialise the map we will use for plotting. The important things to choose here are:1. the projection you wish to use (this may depend on your region of interest, particularly if you are looking at polar data, more information about the different projects is available here: https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html)2. The limits of your map (by default, this will be set to your data limits)3. The resolution of the map coastline.See more information about the options for Cartopy here: https://scitools.org.uk/cartopy/docs/latest/(If you don't have cartopy installed, you can type "conda install -c scitools/label/archive cartopy" in your command prompt, with the anaconda distribution).We start by defining a figure (line 1 below) and then defining a map projection (line 2). All mapping instructions are taken care of using our map object 'm'. Now we make the plot (this may take some time to draw!) You can also save the figure using the code below (this will save in the folder where you are running the code, if you want to save elsewhere you need to specify the path).<jupyter_code>fig2 = plt.figure(figsize=(20, 20), dpi=300)
m = plt.axes(projection=ccrs.PlateCarree(central_longitude=0.0))
f1 = plt.pcolormesh(LON_subset,LAT_subset,np.ma.masked_invalid(CHL_subset), shading='flat', vmin=np.log10(0.01), vmax=np.log10(50), cmap=plt.cm.viridis)  
m.coastlines(resolution=land_resolution, color='black', linewidth=1)
m.add_feature(land_poly)
g1 = m.gridlines(draw_labels = True)
g1.xlabels_top = False
g1.xlabel_style = {'size': 16, 'color': 'gray'}
g1.ylabel_style = {'size': 16, 'color': 'gray'}
cbar = plt.colorbar(f1, orientation="horizontal", fraction=0.05, pad=0.07, ticks=[np.log10(0.01), np.log10(0.1),np.log10(0.5), np.log10(1),np.log10(3),np.log10(10),np.log10(50)]) 
cbar.ax.set_xticklabels(['0.01','0.1','0.5','1','3','10','50'], fontsize=20) 
cbar.set_label('Chlorophyll, mg.m$^{-3}$', fontsize=20)
plt.title('OLCI [CHL_NN] mg.m$^{-3}$', fontsize=20);
plt.show()
fig2.savefig('OLCI_CHL_spatial_demo_no_flags.png', bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>However, this data is not flag masked. This means that we may have data that is subject to glint, or cloud, or a variety of other conditions that variously undermine quality. So, lets apply some flags. We are going to flag extensively, removing all data that corresponds to the following conditions...<jupyter_code>flags_we_want = ['CLOUD', 'CLOUD_AMBIGUOUS', 'CLOUD_MARGIN', 'INVALID', 'COSMETIC', 'SATURATED', 'SUSPECT',
                 'HISOLZEN', 'HIGHGLINT', 'SNOW_ICE', 'AC_FAIL', 'WHITECAPS', 'ANNOT_ABSO_D', 'ANNOT_MIXR1',
                 'ANNOT_TAU06']

file_name_flags = 'wqsf.nc'
FLAG_file = xr.open_dataset(os.path.join(input_root,input_path,file_name_flags))
# get all the flag names
flag_names = FLAG_file['WQSF'].flag_meanings.split(' ')
# get all the flag bit values
flag_vals = FLAG_file['WQSF'].flag_masks
# get the flag field itself
FLAGS = FLAG_file.variables['WQSF'].data
FLAG_file.close()

# make the flag mask using the function we defined above "flag_data_fast"
flag_mask = flag_data_fast(flags_we_want, flag_names, flag_vals, FLAGS, flag_type='WQSF')
flag_mask = flag_mask.astype(float)
flag_mask[flag_mask == 0.0] = np.nan

# subset the flag mask
FLAG_subset = flag_mask[row1:row2:grid_factor, col1:col2:grid_factor]<jupyter_output><empty_output><jupyter_text>And now we apply the flag data to our data and plot again...<jupyter_code>CHL_subset[np.isfinite(FLAG_subset)] = np.nan

fig3 = plt.figure(figsize=(20, 20), dpi=300)
m = plt.axes(projection=ccrs.PlateCarree(central_longitude=0.0))
f1 = plt.pcolormesh(LON_subset,LAT_subset,np.ma.masked_invalid(CHL_subset), shading='flat', vmin=np.log10(0.01), vmax=np.log10(50), cmap=plt.cm.viridis)  
m.coastlines(resolution=land_resolution, color='black', linewidth=1)
m.add_feature(land_poly)
g1 = m.gridlines(draw_labels = True)
g1.xlabels_top = False
g1.xlabel_style = {'size': 16, 'color': 'gray'}
g1.ylabel_style = {'size': 16, 'color': 'gray'}
cbar = plt.colorbar(f1, orientation="horizontal", fraction=0.05, pad=0.07, ticks=[np.log10(0.01), np.log10(0.1),np.log10(0.5), np.log10(1),np.log10(3),np.log10(10),np.log10(50)]) 
cbar.ax.set_xticklabels(['0.01','0.1','0.5','1','3','10','50'], fontsize=20) 
cbar.set_label('Chlorophyll, mg.m$^{-3}$', fontsize=20)
plt.title('OLCI [CHL_NN] mg.m$^{-3}$', fontsize=20);
plt.show()
fig3.savefig('OLCI_CHL_spatial_demo_flags.png', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.089157952669235, 'fraction_numerical': 0.014749587231700605, 'mean_word_length': 4.689104571070757, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 88, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '93809', 'license': 'MIT', 'n_tokens_mistral': 5811, 'n_tokens_neox': 5628, 'n_words': 2110, 'path': 'ocean/Sentinel3/12_OLCI_spatial_interrogation.ipynb', 'repo_name': 'stewartchrisecmwf/wekeo-jupyter-lab', 'size': '18186'}