Updated description of horizontal mapping in operator.

+------------------------------+----------------------------+---------------+-----------------------------------------+---------------------------------+---------------------+

| shape                        | simulation                                                                           | gradient                                              |

+                              +----------------------------+---------------+-----------------------------------------+---------------------------------+---------------------+

|                              | entitie                    | varname       | description                             | entitie                         | varname             |

|                              | entity                     | varname       | description                             | entity                          | varname             |

+==============================+============================+===============+=========================================+=================================+=====================+

| :math:`(n_x,n_y,n_z)`        | :math:`\bar{\mathbf{c}}`   |               | model concentrations                    | :math:`\bar{\mathbf{g}}`        |                     |

+------------------------------+----------------------------+---------------+-----------------------------------------+---------------------------------+---------------------+

The output of the domain-decomposition should be the same as for the serial code.

.. _obsoper-rcfile:

Settings for test program

=========================

The result is computed as a weighted average over grid cells:

.. math::

    \mathbf{G}_p\ \bar{\mathbf{c}}(:,:) ~=~ \left(\ \sum\limits_{k\in K_p} \bar{c}(i_k,j_k)\ w_{p,k}\ \right)/w_p

where the weights represent area in m\ :sup:`2`:

    \mathbf{G}_p\ \bar{\mathbf{c}}(:,:) ~=~ \sum\limits_{k\in K_p} \bar{c}(i_k,j_k)\ w_{p,k}/w_p

* :math:`w_p` is the footprint area of pixel :math:`p`;

* :math:`w_{p,k}` is the fraction of this area that overlaps with the cell :math:`(i_k,j_k)`.

where:

The set :math:`K_p` holds the factioning of pixel :math:`p`.

* set :math:`K_p` includes the indices :math:`k` of the grid cells :math:`(i_k,j_k)`

  that are (partly) overlapped by pixel :math:`p`;

* :math:`w_{p,k}` is the fraction of the pixel area [m\ :sup:`2`] 

  that overlaps with the cell :math:`(i_k,j_k)`;

* :math:`w_p` is the footprint area of pixel :math:`p` [m\ :sup:`2`].

The figure below illustrates how the weights :math:`w_k` are determined.

A footprint is devided into triangles, and the centroids are collected as sample points.

The figure below illustrates how the areas :math:`w_p` and area fractions :math:`w_{p,k}` are estimated.

A footprint is devided into triangles, for which the area is computed using

`Heron's formula <https://en.wikipedia.org/wiki/Heron%27s_formula>`_.

The centroids of the triangless are collected as sampling points.

For each centroid, an originating model grid cell is chosen, in the tutorial code this is simply 

the model grid cell holding the centroid.

The weight of each grid cell in the footprint average is then defined as

the fraction of the area covered by the triangles for the centroids has been assigned to the grid cell.

The area of a pixel that covers a grid cell is then set to the sum of the triangle areas

for which the centroids have been assigned to the grid cell.

The division in triangles is done recursively up to level defined by the user.

Level 0 means no division, and the only centroid is the center of the footprint.

Level 1 means as single division into triangles, leading to 4 centroids.

With level 2, the first 4 triangles are devided into 4 new triangles, leading to 16 centroids in total;

similar, level 3 will lead to 64 centroids.

With level 2, each of the first 4 triangles are devided into 2 new triangles, 

leading to 8 centroids in total;

similar, levels 3, 4, and 5 will lead to 16, 32, and 64 centroids.

The level number that is needed depends on ratio between footprint size and grid cell size;

The level number that is needed depends on the ratio between footprint size and grid cell size;

if the footprints are rather large compared to the model, then a higher level is needed.

For the tutorial code, the level is defined in the settings with::

  ! recursion level for division of footprints into triangles:

  tutorial.mapping.levels            :  3

  tutorial.mapping.levels            :  5

This weights have been used to create the maps on the right in the figure below.

The colors show the sum of the weights that are assigned to a cell during the mapping of the pixels.

At the center of the orbit, this sum could reach values up to 6, which indicates that a grid cell

is covered by about 6 rather small pixels.

For level 0 the sums are rather noisy, but for higher recursion levels the pattern becomes more smooth.

The figure below illustrates the different levels of triangulization.

The maps show the total footprint area associated with the centroids in a grid cell.

For level 0 the area is rather noisy, since in this case the area of a pixel is assigned to

just one grid cell.

For higher recursion levels the pattern becomes more smooth and almost equal to the grid cell area.

.. figure:: figs/footprint-centroids-mapping.png

   :scale: 90 %

.. figure:: figs/footprints_nc4_levels_and_area.png

   :align: center

   :alt: centroids level 0

   *Illustration of mapping from grid cells to footprints.

   Left column shows distribution of (triangle) in 4-sided footprint

   Top row shows distribution of centroids (triangles) in 4-sided footprint

   for different levels of triangularization.

   Maps on the right show total weight of grid cells in demo configuration

   in mapping to S5p orbit. The small figure shows the centers of S5p pixels,

   which are closer to each other at the edges of the swath.*

   Maps below shows per grid cell the sum of the triangle areas for which the centroid is 

   in the cell.*

The CSO code contains a class to perform mapping from a regular grid to footprint area.

Initialize the mapping operator using a grid description::