<jupyter_start><jupyter_text>OLCI spectral plotting     Version: 2.0    Date:    10/04/2019    Author:  <NAME> and <NAME> (Plymouth Marine Laboratory)    Credit:  This code was developed for EUMETSAT under contracts for the Copernicus              programme.    License: This code is offered as free-to-use in the public domain, with no warranty. This routine shows examples of how to use python netcdf libraries to plot OLCI L2 spectra, either for a single point or averaged over and area. We recommend you take a look at the OLCI spatial plotting routine prior to using this one, as it provides more information on some of the steps we are going through here. We start by import libraries for all the functions we need<jupyter_code>%matplotlib inline

import xarray as xr
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import os
import fnmatch
import datetime
import logging<jupyter_output><empty_output><jupyter_text>Here we define a function to calculate spherical distance between points so that we can define polygons<jupyter_code># function
def spheric_dist(lat1,lat2,lon1,lon2,mode="global"):

    '''
      function dist=spheric_dist(lat1,lat2,lon1,lon2)
      compute distances for a simple spheric earth

      input:
      lat1 : latitude of first point (matrix or point)
      lon1 : longitude of first point (matrix or point)
      lat2 : latitude of second point (matrix or point)
      lon2 : longitude of second point (matrix or point)

      output:
      dist : distance from first point to second point (matrix)
    '''

    R = 6367442.76
    # Determine proper longitudinal shift.
    l = np.abs(lon2-lon1)
    try:
        l[l >= 180] = 360 - l[l >= 180]
    except:
        pass
    # Convert Decimal degrees to radians.
    deg2rad = np.pi/180
    phi1    = (90-lat1)*deg2rad
    phi2    = (90-lat2)*deg2rad
    theta1  = lon1*deg2rad
    theta2  = lon2*deg2rad
 
    lat1    = lat1*deg2rad
    lat2    = lat2*deg2rad
    l       = l*deg2rad

    if mode=="global":
        # Compute the distances: new
        cos     = (np.sin(phi1)*np.sin(phi2)*np.cos(theta1 - theta2) + 
                   np.cos(phi1)*np.cos(phi2))
        arc     = np.arccos( cos )
        dist    = R*arc
    elif mode=="regional":
        # Compute the distances: 1 old, deprecated ROMS version - unsuitable for global
        dist    = R*np.arcsin(np.sqrt(((np.sin(l)*np.cos(lat2))**2) + (((np.sin(lat2)*np.cos(lat1)) - \
                  (np.sin(lat1)*np.cos(lat2)*np.cos(l)))**2)))
    elif mode=="local":
        #uses approx for now: 
        x = [lon2-lon1] * np.cos(0.5*[lat2+lat1])
        y = lat2-lat1
        dist = R*[x*x+y*y]^0.5
    else:
        print("incorrect mode")

    return dist<jupyter_output><empty_output><jupyter_text>To help to find your data, please complete the MYPATH variable below with the output generated by the /Configuration_Testing/Data_Path_Checker.ipynb Jupyter notebook in the Configuration_Testing folder.<jupyter_code># e.g. MYPATH = os.path.join("C:/","Users","me","Desktop")
MYPATH = "<please insert your path from Data_Path_Checker.ipynb here, removing the quotes and chevrons>"
#-default parameters------------------------------------------------------------
DEFAULT_LOG_PATH    = os.getcwd()
DEFAULT_FILE_FILTER = '*reflectance*.nc'<jupyter_output><empty_output><jupyter_text>And here we enter the programme proper, setting up our preliminaries and logfile.<jupyter_code>#-main-------------------------------------------------------------------------
# preliminary stuff
logfile = os.path.join(DEFAULT_LOG_PATH,"OLCI_spectral_plot_"+datetime.datetime.now().strftime('%Y%m%d_%H%M')+".log")

# we define a verbose flag to control how much info we want to see. It can also be useful to define a debug flag
# for even more information.
verbose = False
    
# set file logger
try:
    if os.path.exists(logfile):
        os.remove(logfile)
    print("logging to: "+logfile)
    logging.basicConfig(filename=logfile,level=logging.DEBUG)
except:
    print("Failed to set logger")<jupyter_output><empty_output><jupyter_text>Next we define the a box to average our spectra over. If we set the lon values equal, and lat values equal, then we take the spectra from the nearest single point.<jupyter_code># makes box average. Will choose nearest point if latmin == latmax and lonmin == lonmax
lonmin = 70.0
lonmax = 70.25
latmin = 14.0
latmax = 14.25

nearest_flag = False
if lonmin == lonmax and latmin == latmax:
    print('Using nearest point')
    nearest_flag = True
else:
    print('Making box average')<jupyter_output><empty_output><jupyter_text>Now we load our latitude and longitude fields so that we can extract the indices for the box (or point) we want to get the spectra for.<jupyter_code># get the paths and coordinates files    
input_root    = os.path.join(MYPATH,'OLCI_test_data')
input_path    = "S3A_OL_2_WFR____20171226T045629_20171226T045929_20171227T105453_0179_026_076_2700_MAR_O_NT_002.SEN3"
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'][:]
GEO_file.close()

dist_i1 = spheric_dist(latmin,LAT,lonmin,LON)
#J is the X-coord
I1,J1   = np.where(dist_i1 == np.nanmin(dist_i1))

if nearest_flag:
    I1 = I1[0]
    J1 = J1[0]
    I2 = I1+1
    J2 = J1+1
else:
    dist_i2 = spheric_dist(latmax,LAT,lonmax,LON)
    I2, J2   = np.where(dist_i2 == np.nanmin(dist_i2))
    I1 = I1[0]
    J1 = J1[0]
    I2 = I2[0]
    J2 = J2[0]
    
# re-arrange coordinates so that we count upwards...
if J2 < J1:
    J1f = J2
    J2f = J1
else:
    J1f = J1
    J2f = J2
    
if I2 < I1:
    I1f = I2
    I2f = I1
else:
    I1f = I1
    I2f = I2<jupyter_output><empty_output><jupyter_text>The values for the wavelengths of each channel are stored in the xml manifest files, so lets read through that line by line and extract them, as well as the bandwidths.<jupyter_code># get wavelengths from xml file
bands              = []
wavelengths        = []
bandwidths         = []    
xml_file = os.path.join(input_root,input_path,'xfdumanifest.xml')
with open(xml_file, 'r') as input_file:
    for line in input_file:
        if "<sentinel3:band name=" in line:
            bands.append(line.replace('<sentinel3:band name="',"").replace('">','').replace(' ','').replace('\n',''))
        if "<sentinel3:centralWavelength>" in line:
            wavelengths.append(float(line.replace('<sentinel3:centralWavelength>','').replace('</sentinel3:centralWavelength>','')))
        if "<sentinel3:bandwidth>" in line:
            bandwidths.append(float(line.replace('<sentinel3:bandwidth>','').replace('</sentinel3:bandwidth>','')))<jupyter_output><empty_output><jupyter_text>Now lets loop through all of the reflectances, one-by-one and pull out the value over the boox or point for each field, storing them in a list.<jupyter_code># -get the files by band name-------------------------------------------------------------
nc_files=[]
for root, _, filenames in os.walk(os.path.join(input_root,input_path)):
    for filename in fnmatch.filter(filenames, DEFAULT_FILE_FILTER):
        nc_files.append(os.path.join(root, filename))
        if verbose:
            print('Found: '+filename)
        logging.info('Found: '+os.path.join(root, filename))

# get the reflectances
reflectances       = []
reflectance_errors = []
variability        = []
for nc_file in sorted(nc_files):
    varname  = os.path.basename(nc_file).split('.')[0]
    nc_fid   = xr.open_dataset(nc_file)
    reflectances.append(np.nanmean(nc_fid.variables[varname][I1f:I2f,J1f:J2f]))
    variability.append(np.nanstd(nc_fid.variables[varname][I1f:I2f,J1f:J2f]))
    reflectance_errors.append(np.nanmean(nc_fid.variables[varname+'_err'][I1f:I2f,J1f:J2f]))
    nc_fid.close()

tops = [x + y for x, y in zip(reflectances, variability)]
bottoms = [x - y for x, y in zip(reflectances, variability)]<jupyter_output><empty_output><jupyter_text>Finally, we plot the spectra and save.<jupyter_code># print the spectra
fig1 = plt.figure(figsize=(10, 10), dpi=300)
plt.errorbar(wavelengths, reflectances, xerr=bandwidths, yerr=reflectance_errors,color='k',linewidth=1.0)
plt.plot(wavelengths,tops,'r--',linewidth=1.0)
plt.plot(wavelengths,bottoms,'r--',linewidth=1.0)
plt.fill_between(wavelengths,tops,bottoms,color='r',alpha=0.25)
plt.xlabel('Wavelength [nm]', fontsize=16)
plt.ylabel('Water leaving reflectance [sr$^{-1}$]', fontsize=16)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.title('OLCI L2 reflectances', fontsize=16);
fig1.text(.5, -0.05, r'OLCI L2 water leaving reflectance spectra. Black trace shows the mean spectrum,'+'\n'+ \
          r'with +/- 1 standard deviation covered by the red shaded area. Horizontal and vertical'+ '\n' + \
          r'bars show the respective band widths and reflectance errors for each wavelength.' +'\n' + \
          r'Note that water leaving reflectance ($\rho_w$) can be converted to remote sensing' +'\n' + \
          r'reflectance ($R_{rs}$) by dividing by $\pi$.', ha='center')
plt.show()
fig1.savefig('OLCI_spectra_demo.png',bbox_inches='tight');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13698481561822126, 'fraction_numerical': 0.030585683297180043, 'mean_word_length': 3.4225419664268584, 'pattern_counts': {'":': 3, '<': 50, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 5}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '85109', 'license': 'MIT', 'n_tokens_mistral': 3250, 'n_tokens_neox': 3087, 'n_words': 830, 'path': 'OLCI/13_OLCI_spectral_interrogation.ipynb', 'repo_name': 'africamarine-eo/training', 'size': '9236'}