Using output files - ORAC-CC/orac GitHub Wiki

Python module

The Anaconda package installs a Python module called pyorac which contains the Swath class, designed for managing ORAC output files. (If working from source, simply run all code below from the tools directory or add that directory to your PYTHONPATH variable.) To use initialise the class,

from pyorac.swath import Swath

orac_obj = Swath("/path/to/file.nc")
with Swath("/some/other/file.nc") as orac_two:
    # Some code

If the primary and secondary output files are in the same directory, the class will automatically open both, allowing you to access all variables from a single object. Variables contents can accessed as if the object was a dictionary, e.g. data = orac_obj["cot"] will store the contents of the field cot as a numpy array in data. Access to the variable's attributes would be done through orac_obj.get_variable("cot").

A number of derived variables are available through class properties:

• Swath.cldmask returns the binary cloud mask;
• Swath.cldflag returns the cloud clearing field (see below);
• Swath.ang returns the 550-870 Angstrom exponent and Swath.ang_unc it's uncertainty;
• Swath.rs return the surface reflectance (a 3D array, as this has a spectral component);
• Swath.qcflag returns the quality control flag.

The following methods are available:

• Swath.flag_map() returns a dictionary mapping the values of a bitmask to strings describing their meaning;
• `Swath.mask_cloud()' returns a Boolean array that is True where the aerosol cloud flag is non-zero, suitable for extracting aerosol pixels from a field;
• Swath.mask_clear() returns a Boolean array that is True where the cloud mask is non-unity, suitable for extracting cloud pixels from a field;
• Swath.map() is a wrapper for matplotlib.pcolormesh to plot the satellite swath on a cartopy axes. The first argument is the axes to plot on and the second is the data field to plot, with all keywords passed directly to pcolormesh;
• Swath.false() is a wrapper for matplotlib.pcolormesh to plot false-colour images on a cartopy axes. The first argument is the axes to plot on and the second is the data field to plot (for which the first dimension is used to derive RGB, respectively, though the order keyword may be used to specify an alternative subset).

Cloud flagging

There are multiple cloud masks available within ORAC. Roughly, cloud masks indicate pixels with a successful cloud retrieval (i.e. a conservative cloud flag) while cloud flags indicate pixels that may contain cloud (i.e. a conservative clear-sky flag). Swath.cldmask is the official Cloud Mask of the Cloud CCI program, generated by a neural network. Swath.cldflag is an implementation of the cloud clearing algorithm used by ORAC within the Aerosol CCI program.

Other classes

pyorac.swath also contains the Flux and MergedSwath classes. The former manages the outputs of the broadband flux code while the latter allows the simultaneous opening of 1km cloud and 10km aerosol data (when a collocation file is available). These are not yet documented thoroughly and are under development. The plotting functionality is inherited from pyorac.mappable.Mappable if the user has an interest in using those routines independently.

Examples

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
from pyorac.swath import Swath

def primary_map(ax, values, cb_pos, label, perc=(5,95)):
    """Make a nicely formatted map of a variable."""
    # Map values
    vmin, vmax = np.percentile(values.compressed(), perc)
    im  = orac.map(ax, values, rasterized=True, vmin=vmin, vmax=vmax)
    # Formatting appearance
    ax.set_aspect("auto")
    ax.coastlines("10m")
    gl = ax.gridlines()
    # Colourbar
    cb = fig.colorbar(im, ax=ax, label=label)
    return cb

f = "/data/povey/data/norm/MEXICO/ESACCI-L2-CLOUD-CLD-AATSR_ORAC_Envisat_200806031603_R1787.primary.nc"
with Swath(f) as orac:
    fig, ax = plt.subplots(1, 1, subplot_kw=dict(projection=ccrs.PlateCarree()))
    _ = primary_map(ax, orac["aot550"], [.05, .71, .02, .17],
                        "Aerosol optical depth")
    plt.show()

IDL Script

The ORAC outputs are in NCDF format, for which I/O functions are included in the standard IDL release (see http://www.exelisvis.com/docs/NCDF_Overview.html for more information). The repository contains IDL routines for producing evaluation plots for debugging purposes. One of these routines, http://proj.badc.rl.ac.uk/orac/browser/trunk/idl/ncdf_obtain.pro, is generally useful for reading a data field of known name (the field names are given in here) as it applies the scale factor, offset, and missing value attributes.

The following is an example of reading the cloud optical thickness from a file (the variable `filename' should be set previous to this code) and making a basic plot using routines found in the repository.

   ;; Open file
   fid = NCDF_OPEN(filename)

   ;; Read out desired data fields
   lat = NCDF_OBTAIN(fid, 'lat')
   lon = NCDF_OBTAIN(fid, 'lon')
   cot = NCDF_OBTAIN(fid, 'cot')

   ;; Close file
   NCDF_CLOSE, fid

   ;; Plot the data on a map projection
   MAPPOINTS, cot, lat, lon

That makes use of two EODG routines. MAPPOINTS is our procedure to plot a 2D array as points over a map and can be distributed on request (it's long and no longer actively supported). NCDF_OBTAIN is,

;+
; NAME:
;   NCDF_OBTAIN
;
; PURPOSE:
;   Fetch a field from an open NetCDF file, applying the appropriate scale
;   factor, offset, and fill value.
;
; CATEGORY:
;   ORAC plotting tools
;
; CALLING SEQUENCE:
;   field = NCDF_OBTAIN(file_id, name [, fill])
;
; INPUTS:
;   file_id = An ID number returned by NCDF_OPEN.
;   name    = The name of the field to be opened.
;
; OPTIONAL INPUTS:
;   None.
;
; KEYWORD PARAMETERS:
;   None.
;
; OUTPUTS:
;   field   = The requested data.
;
; OPTIONAL OUTPUTS:
;   fill    = The fill value applied to the data.
;
; RESTRICTIONS:
;   None.
;
; MODIFICATION HISTORY:
;   15 Jul 2014 - ACP: Initial version ([email protected]).
;    2 Feb 2015 - ACP: Change selection of fill value to be based on output format.
;-
FUNCTION NCDF_OBTAIN_FILL, in
   ON_ERROR, 2
   COMPILE_OPT HIDDEN
   ;; determine appropriate fill value
   case SIZE(in,/type) of
      5: fill = !values.d_nan
      4: fill = !values.f_nan
      3: fill = -2147483647l
      2: fill = -32767
      1: fill = 255
      else: MESSAGE, 'Fill value not meaningful for this field. '+ $
                     'Use a different routine.'
   endcase
   return, fill
END

FUNCTION NCDF_OBTAIN, fid, name, fill, sc, off, id=id
   COMPILE_OPT HIDDEN, LOGICAL_PREDICATE, STRICTARR, STRICTARRSUBS

   ;; read data array from NCDF file and apply necessary scalling
   vid = KEYWORD_SET(id) ? name : NCDF_VARID(fid,name)
   if vid lt 0 then RETURN, -1
   NCDF_VARGET,fid,vid,data
   vq=NCDF_VARINQ(fid,vid)
   nc_fill=!values.f_nan

   for aid=0,vq.natts-1 do begin
      nm=NCDF_ATTNAME(fid,vid,aid)
      case nm of
         '_FillValue': begin
            NCDF_ATTGET,fid,vid,nm,nc_fill
         end
         'scale_factor': NCDF_ATTGET,fid,vid,nm,sc
         'add_offset': NCDF_ATTGET,fid,vid,nm,off
         else:
      endcase
   endfor

   ;; replace fill value
   p_fill=WHERE(data eq nc_fill,nfill,comp=p_valid,ncomp=nvalid)

   if KEYWORD_SET(sc) then begin
      fill = NCDF_OBTAIN_FILL(sc)
      if KEYWORD_SET(off) then begin
         out = REPLICATE(sc,SIZE(data,/dim))
         if nvalid gt 0 then out[p_valid] = off + sc*data[p_valid]
      endif else begin
         out = REPLICATE(sc,SIZE(data,/dim))
         if nvalid gt 0 then out[p_valid] = sc*data[p_valid]
      endelse
   endif else begin
      if KEYWORD_SET(off) then begin
         fill = NCDF_OBTAIN_FILL(off)
         out = REPLICATE(off,SIZE(data,/dim))
         if nvalid gt 0 then out[p_valid] = off + data[p_valid]
      endif else begin
         fill = NCDF_OBTAIN_FILL(data[0])
         out = data
      endelse
   endelse
   if nfill gt 0 then out[p_fill] = fill

   RETURN, out

END