Hello ARTS

Posted on August 21, 2017
Tags: ARTS, python

To simplify development and debugging of the retrieval code that I am currently working on in ARTS, I started writing a python interface. After a few (more) hours of digging into ARTS and fixing segfaults, the new python interface can be be used as a fully interactive, bi-directional interface to ARTS providing

In this post, I’ll write up a short Hello World to demonstrate the capabilities of the interface so far.

Setup

The interface is available from the development version of typhon and requires the ARTS C API, which has been added to the development version of ARTS. Building the C API must be enabled using:

cmake ... -DENABLE_C_API=1

The python interface needs access to the arts_api.so shared library, which is located in the src subfolder of the ARTS build tree. The interface expects the location of the ARTS build tree to be provided in the ARTS_BUILD_PATH environment variable.

The Workspace Class

The main functionality of the interface to ARTS is provided by the Workspace class contained in the the workspace sub-package of typhon. A Workspace object represents an ongoing ARTS simulation:

from typhon.arts.workspace import Workspace

ws = Workspace()

Executing Controlfiles

Similar to an INCLUDE statement in an ordinary controlfile, external controlfiles can be executed on a workspace object to perform common setup operations:

ws.execute_controlfile("general.arts")
ws.execute_controlfile("continua.arts")
ws.execute_controlfile("agendas.arts")
ws.execute_controlfile("planet_earth.arts")

Currently the python interface recursively searches for controlfiles with the given name in the current directory and the arts include path given by the ARTS_INCLUDE_PATH environment variable.

In order for the ARTS runtime to be able to resolve INCLUDE statements in an executed controlfile, the interface sets the arts search path to the controlfiles subfolder of the path given in ARTS_INCLUDE_PATH. The data path of the runtime is set by the interface to the controlfiles/testdata folder.

Both, the search path and the data path of the ARTS runtime can also be set manually using the include_path_add and data_path_add methods provided in typhon.arts.workspace.api.

Executing Workspace Methods

ARTS Workspace methods are accessible as class methods of each Workspace object. In IPython also autocompletion should work. Calls to super-generic WSMs are resolved by the interface:

ws.Copy(ws.abs_xsec_agenda,         ws.abs_xsec_agenda__noCIA)
ws.Copy(ws.propmat_clearsky_agenda, ws.propmat_clearsky_agenda__OnTheFly)
ws.Copy(ws.iy_main_agenda,          ws.iy_main_agenda__Emission)
ws.Copy(ws.iy_space_agenda,         ws.iy_space_agenda__CosmicBackground)
ws.Copy(ws.iy_surface_agenda,       ws.iy_surface_agenda__UseSurfaceRtprop)
ws.Copy(ws.ppath_agenda,            ws.ppath_agenda__FollowSensorLosPath)
ws.Copy(ws.ppath_step_agenda,       ws.ppath_step_agenda__GeometricPath)

ws.AtmosphereSet1D()
ws.IndexSet(ws.stokes_dim, 1)

Arguments to workspace methods can be given directly as python objects. This works for the following combinations of ARTS groups and python types:

ARTS Python
Index int
Numeric numpy.float64
String str
ArrayOfString [str]
Vector numpy.ndarray
Matrix numpy.ndarray

If positional arguments are used to call a WSM, at least all generic outputs and non-generic inputs must be given in order. In the call to VectorNLinSpace below, this is used to overwrite the nelem input with a value provided directly from python. This avoids having to set nelem first using the IndexSet workspace method.

# Frequency and pressure grids
f_start = 110.436e9
f_end   = 111.236e9
n_f     = 801
n_p     = 91

ws.VectorNLinSpace(ws.f_grid, n_f, f_start, f_end )
ws.VectorNLogSpace(ws.p_grid, 361, 500e2, 0.1 )

# Spectroscopy
ws.abs_speciesSet(species=[ "O3" ])
ws.abs_lineshapeDefine("Voigt_Kuntz6", "VVH", 750e9)
ws.abs_linesReadFromArts("ozone_line.xml", 0.0, 999e9)
ws.abs_lines_per_speciesCreateFromLines()

If not explicitly given, non-generic inputs are taken from the method definition in methods.cc. Generic inputs can be provided using named arguments with the same names as in the definition in methods.cc. Also default parameters work:

ws.AtmRawRead( basename = "tropical" )
ws.AtmFieldsCalc()

In IPython, documentation of workspace methods can be accessed using the built-in documentation:

$ ws.AtmFieldsCalc()?

Signature: ws.AtmFieldsCalc(*args, **kwargs)
Docstring:
Interpolation of raw atmospheric T, z, VMR, and NLTE T fields to
calculation grids.

An atmospheric scenario includes the following data for each
position (pressure, latitude, longitude) in the atmosphere:
   1. temperature field
   2. the corresponding altitude field
   3. vmr fields for the gaseous species
...

Workspace Variables

Symbolic representations of all workspace variables are provided in the typhon.arts.workspace.variables module. If one doesn’t mind messing up the global namespace, these can be used to write variable names directly in function calls:

from typhon.arts.workspace.variables import *

ws.MatrixSet(z_surface, 10e3 * np.ones((1,1)))
ws.VectorSet(lat_true, np.array([10.0]))
ws.VectorSet(lon_true, np.array([123.0]))

Workspace variable can also be accessed as attributes of a workspace object. This will dynamically look up this attribute in the current workspace and return a symbolic representation of the variable.

ws.lat_true

If a workspace variable is associated to a workspace, it’s value can be accessed using the value() workspace method. This will return a python object representing the value of the workspace variable.

$ ws.lat_true.value()
  array([ 10.])

The print() method of a workspace variable prints out the value of the variable in the given workspace using the ARTS Print(…) WSM:

$ ws.lat_true.print()
  10

If the WSV is of type Vector, Matrix or Tensor it will implement the numpy array interface. This means this variable can be used just like a numpy.ndarray (using numpy.asarray) without copying the data.

import matplotlib.pyplot as plt
import numpy as np

ws.atmfields_checkedCalc()
ws.NumericSet( p_hse, 100e2 )
ws.NumericSet( z_hse_accuracy, 0.5 )
ws.z_fieldFromHSE()

plt.plot(ws.p_grid, np.asarray(ws.z_field).ravel())
plt.ylabel(r"z $[m]$")
plt.xlabel(r"p $[Pa]$")

Python Scripts as Controlfiles

As you might have already guessed, this Hello World is actually the transcript of an ARTS controlfile that I use for testing the OEM implementation. So let’s move forward towards a simple radiative transfer simulation:

ws.NumericSet( ppath_lmax, -1.0 )
ws.StringSet( iy_unit, "RJBT" )
ws.MatrixSet( sensor_pos, 15e3 * np.ones((1,1)))
ws.MatrixSet( sensor_los, 60 * np.ones((1,1)))
ws.VectorSet( sensor_time, np.zeros((1)))

# Deactive parts not used and perform all remaining tests
#
ws.cloudboxOff()
ws.sensorOff()
ws.abs_xsec_agenda_checkedCalc()
ws.propmat_clearsky_agenda_checkedCalc()
ws.atmgeom_checkedCalc()
ws.cloudbox_checkedCalc()
ws.sensor_checkedCalc()
ws.jacobianOff()
ws.yCalc()

Manipulating Workspace Variables

The numpy array types returned by calling numpy.asarray on a workspace variable that is a vector, matrix or tensor effectively returns a view on this variable and can thus be manipulated directly using any operation that can be used on numpy.ndarrays:

plt.figure()

vmr = np.asarray(ws.vmr_field) # View on vmr field.
d_vmr = 0.5 * vmr

plt.plot(ws.f_grid, ws.y, label=r"$vmr = vmr_0$")
for i in range(4):
    vmr += d_vmr # Add perturbation to ARTS vmr field.
    ws.yCalc()
    plt.plot(ws.f_grid, ws.y, label=r"$vmr = vmr_0 + " + str(f + 1) + "\Delta vmr$")

plt.xlabel("f [Hz]")
plt.ylabel("y [RJBT]")
plt.legend()

The above snippet, for example, performs the simulation of the ozone line at 110 GHz over a range of different volume mixing ratios and combines the result in a single plot:

Conclusions and Disclaimer

So far, I have developed this interface solely with the goal of running simple retrievals and had not had the time to test it on many other tasks. Hence there are probably a million ways to break the interface, which will most likely end up in a segfault sooner or later.

However, I thought that this may be useful for some users of typhon and maybe even people who want to learn how to use ARTS. Furthermore, this interface opens up quite a few interesting possibilities of using ARTS programmatically, for example to test different retrieval methods or integration with statistical models.