Data Calculator Instructions

Content-

• Basic Information

• Conceptual Approach
• Anchor and Chain Datasets
• Entering an dataset symbols (variables)
• Operators and basic functions

• Datasets with Time Series Layers

• Single Layer Referencing
• Range or Multiple of TS Layer Referencing
• Time Series Temporal Resolutions

• Reduction of Dimensions in the Output
• Map Styling and Other Options
• Custom Statistical and Other Functions
• Examples (the Most Useful Section!)
• Appendices

• Country FIPS codes
• US State Codes
• River Basin codes
• Continent codes

Basic Information

Data Calculator allows a user to execute algebraic and logical expressions/equations over any number of the gridded datasets in the system on the fly and display results in a form of maps, graphs, and just numbers with an option to download these results as data files in NetCDF (gridded datasets) or ASCII CSV (spreadsheet) format. The calculations are performed on the grid cell level over the full grid extent of the anchor dataset (see definition below).

"Equation" input form-
The equation has to be entered in this  input form and it can contain basic and advanced algebra operators, common mathematical and statistical functions, conditional operators, etc. These are explained below in the Operators section. A set of special instructions could be used to control the calculation results and graphical design of the result visualization by the map and the histogram. These instructions are accessible via form input options.

Conceptual Approach

The calculation equation submitted to the server gets parsed, valideated and then evaluated using the powerful Perl PDL (Perl Data Language) scripting language and data management framework of RIMS (Rapid Integrated Mapping System) developed by WSAG UNH. PDL operates over multi-dimensional arrays of binary data with multi-threading over assigned dimensions and is computationally very efficient and commonly outperforms single threaded C or FORTRAN compiled codes. The datasets entered in the equation are read from the source data files, re-gridded to anchor dataset extents and re-packed to the following dimensions order-

(x, y, t)

where the last dimension t is time and which always has size 1 for non time series datasets (e.g. elevation) or N for the number of time series layers. Accordingly t > 1 for multi-layer time series datasets. First two dimensions x and y represent a two-dimensional data layer.

Understanding primitive dimensions of the datasets in the equation is important for using operators and functions that change (reduce) dimensionality of the return/output data arrays.

Anchor and Chain Datasets

Let us define anchor and chain datasets

• The "anchor" dataset is the very first one used in the equation. And, thus, there is only one anchor dataset in the equation.
• The "chain" datasets are all others, if any.

Anchor dataset controls and determines the following grid properties for all chain datasets in the equation-

• Layer dimensions (x, y)
• Layer geographical projection, extent, and grid cell resolution
• Mask of nodata grid cells

The anchor dataset is taken "as is" for the source data file(s) while the chain datasets are re-projected, re-gridded, and nodata re-masked to match the anchor using GDAL library and tools. So, if your anchor dataset has, say, only continental US spatial domain extents and the slave dataset has originally global extent, the latter will be reduced to continental US extent and the output (if grid) will be in it as well.

After regridding and copying of nodata mask from anchor to a chain dataset, the latter may still have nodata pixels in the “good” data mask of the anchor. These are patched by a valid values as shown in the figure below-

RIMS has three methods for patching the nodata pixels in the chain datasets:

1. Fixed value (common) - all nodata mismatch pixels are patched with the same value, e.g. zero.
2. Secondary dataset (least common) - pixels are patched by values from the alternative secondary dataset (grid), e.g. high quality climate drivers of continental US domain are patched with lower quality Global climate driver layers for, say, North America extents of the Anchor dataset spatial domain. This nodata patching method is not available for the Data Calculator.
3. Nodata erosion (most common) - The nodata pixel adjacent to valid pixels are patched by their average so all nodata pixels around the boundary of valid pixels are patched at a time. The patching process continues until all mismatched pixels around the boundary of valid pixels are done. The “islands” of the mismatched pixels inside continuous areas of nodata in the chain dataset (e.g. the pixel in the upper right corner of the illustration above) are patched by the global average of all valid “continental” pixels.

Tip to clip spatial domain by an anchor dataset:

 An anchor dataset can be used to clip calculation spatial extent to match the domain of the anchor. This can be done by applying a dummy statement at the beginning of the calculation equation. E.g. -

Runoff_d{2000-01-01}->isgood * ( lin_reg_slope(M2_t2m_y{}) )

The Runoff_d dataset has an extent of the Ob river basin only while the M2_t2m_y dataset (MERRA2 annual temperatures at 2 m) has a Global extent. Using Runoff_d as anchor results in clipping the calculations to the Ob river area as illustrated below -

Entering an dataset symbols (variables)

Equation has to present a valid mathematical and/or logical expression that contains a combination of variables and operators. In the case of this Calculator, variables are numerical gridded datasets as the calculations are applied to each of their grid cells/pixels. Dataset notations (symbols) are linked to the actual datasets that are available for Web mapping visualization and hosted on WSAG servers and NAS systems. The dataset symbols that can be used in the equation input are subdivided into three categories-

• Server side pre-set dataset symbols grouped by a theme (e.g. Land, Climate, etc.)
• User defined dataset symbols on the client side that are applied to the primary dataset currently displayed in the map view area. These can be accumulated as a user switches from one map to another and keeps assigning symbols to each or the ones needed for the calculations.
• Generic “all-the-time” symbol “mData” that is linked to the presently viewed dataset.

The server preset dataset symbols are grouped by a theme and listed in the table on the left of the Data Calculator area. See table circled in red ink below-

Various symbol themes/groups are accessible from the radio button switches as shown above in the top row. For example, the Land theme is circled in red ink which is displayed in the symbol table. For the user's convenience, mouse clicking on the symbols in the Dataset Table automatically pastes them to the active Equation input at the active caret position (see Elev symbol circled in grey ink that is pasted in with the link click) . This ensures that the symbols are not misspelled as these are also case sensitive. Of course, a user can type in dataset symbols directly. Also, clicking on the Input dataset name makes it to be displayed on the map above (see “Elevation” circled in the illustration above).

User defined dataset symbols for the datasets currently displayed on the map can be added by clicking “define” link next to the default dataset symbol-

After it is typed in, clicking the “Add” button validates it and adds to the list of User Defined symbols that can be displayed in the symbol list table described above and accessible through clicking “User Defined” radio button (see illustration above).

And finally, the generic “all-the-time” symbol “mData” is available as well which is always available and clicking on it pastes it to the equation at the current caret position.

Note, if a pre- or user-defined symbol for the dataset currently displayed on the map is always displayed as shown here-

For example, a very simple calculation to convert elevation from meter to kilometer units could be done by this equation Elev/1000 where symbol "Elev" represents the elevation dataset.

Operators and Basic Functions

Common algebraic and logical operators are those used by Perl PDL which are the same as in any most higher level programming languages. Few of these have PDL specific exceptions. In addition all default and some custom installed PDL function libraries are available as well. See tables below-

Basic Algebraic Operators-

Logical Operators-

* Note- Logical “and” and ”or” operations in the Calculator equations must be bitwise operations

Logical Functions specific to PDL-

Notes-

• PDL does not have ternary operator, so RIMS has a simple functions for these
• Use condition_slice() function in cases where NaN can result in b or c inputs, e.g. division by zero: condition_slice(PopDensity != 0, LandCost/PopDensity, 0)

Some Common Mathematical Functions-

Some Common Array Operators/Functions-

Datasets with Time Series Layers

Single Layer Referencing

Datasets with Time Series (TS) layers are indicated by a star (*) symbol in the Data Table. Each layer can be accessed via TS layer identifier in a form {YYYY-MM-DD}. For example, a difference in Sea Surface Temperature between two dates (summer and winter) could be calculated by this equation-

SeaST{2007-07-15} - SeaST{2007-12-15}

Note- In the case of this equation a user can enhance the resulting map color range by setting  Map Styling options for Min = -10 and Map Options -> Max = 20.

If no TS layer identifier is provided for a given time series dataset than the latest (last) data layer is used in the equation, e.g. equation

SeaST * 9/5 + 32

converts sea surface temperature from Celsius to Fahrenheit for the latest data layer (near real-time) in sea surface temperature (SeaST) dataset. And if a user want to do the same for some specific date, say November 12, 2007 then it could be done with this formula-

SeaST{2007-11-12} * 9/5 + 32

Range or Multiple of TS Layer Referencing

A syntax to call a set or range of Time Series (TS) layers is provided by a comma and/or double dot operator (range) in the form of

SeaST{YYYY-MM-DD, YYYY-MM-DD..YYYY-MM-DD, etc.}

where, in the case of a double dot range operator, the first data has to be earlier than the second date. Note that TS layers are added to the last dimension of the SeaST data array-

(x, y, N)

where N is a number of TS layers from a user entered set or range operators.

The shortcut or wrapper for the range operator {} is also used in the equation input form as shown in the illustration below which is equivalent to a longer way of entering the equation-

SeaST{2020-07-01 .. 2020-07-31}->mv(2,0)->average

Note, this equation returns the average monthly SeaST temperature for July, 2020.

The double dot operator returns an array of layers so that array functions could be used on them (such as min(), max(), average(), sum(), map(), etc.). For example, a number of days with surface precipitation over 100 mm/day in the month of July, 2020 could be calculated with this equation-

(PrecipN{2020-07-01..2020-07-31} > 100)->mv(2,0)->sumover

Time Series Temporal Resolutions

Different time series temporal resolutions are referenced by their specific form of referencing which are indicated in the table below-

The logic of the date syntax is that it is the same format for all (XXXX-XX-XX) where zeroes denote something that is not available as, say, days in the monthly datasets.

Reduction of Dimensions in the Output

As it was noted in the conceptual approach of the Data Calculator, the input data has always three dimensions

(x, y, t)

where the last (time) dimension can have value 1 for for single layer datasets (e.g. elevation). Some calculations may not reduce any dimensions of the output, e.g. unit conversion. But many types of calculations commonly collapse the third (time) dimension to value 1 (one layer) or completely eliminate it, e.g. calculation of min, max, average, etc. values of a time series data. In some calculations x or y dimensions can be collapsed as well which are those where output is a longitudinal or latitudinal cross-sections of the data, e.g. Global percent of land by the latitude. All aforementioned cases are allowed and they result in the “grid” “table” or “number” types of output. See table below.

Table notes:

•  t notation indicates t>1 (multi-layer), and it is replaced by 1 for single-layer datasets.
• TS stands for Time Series which is for all cases of t>1.
• v is always v>1, otherwise the output type is Number.

Calculation output results are displayed back to the user depending of the output type (indicated in the table above) as follows-

• Grid - Map view with full set of features and controls.
• Table - Graph window popup with select form for each v or t value.
• Number - A popup window with output number displayed.

If a popup window is closed or the user needs to access metadata for the calculations, download calculation results data (file) or init metadata file for Grid output, these are alway available for the last performed calculation by clicking “Last Calculation Info/Download” button in the lower right corner of the Calculator-

Map Styling and Other Options

Map Styling allows an option to choose color palette, its min/max range, and whether it should be in log10 scale. These have an effect only for the “Grid” output type which is displayed in the map view window. The default color palette is “Rainbow” and the default min max values are taken from the range of calculated values. The choice of color palette is quite useful as most of Precipitation calculations look much better in the “DarkBlue” colors, vegetation or irrigation in the “DarkGreen”, temperatures in “SST”, and “all over” outputs are better in “ROYG” color palette, etc.

Presently, map server processing options include controls only on “noScale” and “noOffsite” inputs, but we are planning on including a full set of MT “Processing” keys in the near future.

Other Options include choices for reading of the Chain datasets (see “Anchor and Chain Datasets” section above). These are-

• Chain dataset resample reading option which can be GDAL interpolation for Nearest Neighbor or Bilinear (default) methods. While the bilinear default is certainly best for most of calculations, but some may need to use the “Nearest Neighbor”, e.g. cases with the use of polygon masks with indices such as countries, states, land cover types, etc. where interpolation (bilinear) between two integer source values is not a wrong choice as it can lead to unwanted results.
• Nodata patch value affects how the chain dataset nodata values are converted to data values to match valid pixels of the anchor dataset. Note that the anchor dataset nodata pixel mask is copied over to all chain datasets first. Depending on the input for the nodata value (see illustration)
          
  two possible patching methods are used:

• Fixed value (default is zero), or
• Nodata erosion (both described in the section Anchor and Chain Datasets), if the “Nodata patch value” input is empty (no number).

• Trim nodata margins option removes nodata columns and rows around valid data of the calculation gridded output. That makes the resulting files smaller and more useful for re-use outside of the Web application. See illustration of that option location in the Calculator-
  

Custom Statistical and Other Functions

These are created for user convenience to perform some common statistical and other calculations that may require a long statement or set of statements in the Equation Input form. For example, calculation of linear regression slope (say, for climate warming trend in C/year) can be done using this input (no line wrap)-

my @dim=era5_t2m_y{}->dims;

my $x=sequence($dim[2])->dummy(0,$dim[1])->dummy(0,$dim[0]);

my $y=era5_t2m_y {};

($dim[2]*($x*$y)->mv(2,0)->sumover - $x->mv(2,0)->sumover*$y->mv(2,0)->sumover)/($dim[2]*($x**2)->mv(2,0)->sumover - ($x->mv(2,0)->sumover)**2)

But it is more convenient to do it using a pre-built custom function that does exactly the same as the lengthy input above-

 lin_reg_slope(era5_t2m_y{})

The list of such custom functions is listed in the Table below-

Table of Custom Functions-

These function are listed and available for direct paste in here-

Please, let me know, if you want to add more custom functions

Examples (the Most Useful Section!)

 As it was noted in the  “Basic Information” section, the Data Calculator is using algebraic and logical operators, common functions and equation syntax of Perl PDL (Perl Data Language) where a number of links to useful in-depth documentation and information can be found.

First, let us parse and explain some fundamentals of a few equations for some most common calculation tasks.

Global average of time series layers, e.g. calculation of average Global temperature (single number) from time series of annual 1980-2019 MERRA2 air temperature at 2 m. This is the equation to do it in variants: (1) useless unweighted simple average, and (2) area area weighted average which is very important to do since the source dataset is in the Geographical projection -

1.  M2_t2m_y{}->avg

2. (M2_t2m_y{}*CellArea)->sum / (CellArea->sum * M2_t2m_y{}->dim(2))

where time series range {} is entered as {1980-00-00..2019-00-00} in the form below the equation as shown in this illustration

Here, for simplicity, we are going to explain and parse the equation (1) where the very first part of the equation M2_t2m_y{} is the input array (x,y,t) with actual dimensions of (576, 361, 40) , i.e. 40 time series layers of annual temperatures and we need to calculate an average of all pixel values in all time series layers. The PDL function ->avg does an average over all dimensions resulting in a single number (see Operators and Basic Functions section above)

Average layer of time series layers, e.g. calculation of long term average temperature (single grid layer) from time series of annual 1980-2019 MERRA2 air temperature at 2 m. This is the equation to do it

M2_t2m_y{}->mv(2,0)->average

where time series range {} is entered as {1980-00-00..2019-00-00} in the form below the equation as shown in the illustration for the previous example. Note, the cell area weighing is not needed here since the output is grid.

Similarly to the first example above, the very first part of the equation M2_t2m_y{} is the input array (x,y,t) with actual dimensions of (576, 361, 40), i.e. 40 time series layers of annual temperatures and in this example we need to calculate an average of 40 values in each pixel over t dimension. The PDL function ->average is the one to use as it does averaging (see Operators and Basic Functions section above), but only over the first dimension (which is x) of the (x,y,t) array. So, we need to reorder the dimensions (#0,#1,#2) so that the time dimension (#2) moves to the place of the first dimension (#0). PDL function for reordering is ->mv(2,0) which is moving dimension #2 to the place if #0. So the result of M2_t2m_y{}->mv(2,0) part of the equation yields an array of (t,x,y). Now, the function ->average makes an average of the first dimension t and reduces dimensionality to (x,y) making an output grid layer that is displayed on the map. The output is saved to NetCDF file which can be downloaded by clicking on the “Last Calculation Info/Download” button.

Global average of each time series layer, e.g. calculation of global average temperature (single number for each layer) from time series of annual 1980-2019 MERRA2 air temperature at 2 m. This is the equation to do it in both (1) unweighted (useless), and (2) weighted forms

1.  M2_t2m_y{}->average->average

2. (M2_t2m_y{}*CellArea)->sumover->sumover / CellArea->sum

where time series range {} is entered as {1980-00-00..2019-00-00} in the form below the equation as shown in the illustration for the previous example.

Similarly to the previous examples above, the very first part of the equation M2_t2m_y{} is the input array (x,y,t) with actual dimensions of (576, 361, 40). As we have already learned the PDL function ->average is making averaging over the first dimension x, so that the resulting array of reduced dimensionality becomes (y,t). And repeating it over (y,t) with ->average makes an average of y dimension, making the result to have only one remaining dimension (t) which is the array of 40 numbers of average global temperature for each year from 1980 to 2019. This is “graph” output type of the Data Calculator that looks like this -

A. Table of some useful examples of calculation equations
        Note 1: use TS-Range input for {} in the equation, e.g. 1980-00-00 .. 2019-00-00
        Note 2: NN - use “Nearest Neighbour” read method for Chain datasets (e.g. State, Country)

Appendices

Country FIPS codes

These are index IDs in the “Country” input dataset (the first one in the “Other” theme) -

This dataset contains enumerated standard FIPS codes for World countries (see table below) as country index IDs. These IDs are useful for many calculations where a single or list of countries has to be masked or processed.

Usage example: Population-2000 of China where 46 is ID value for China-

(Pop00*CellArea)->where(Country == 46)->sum

Table of country FIPS codes and ID values-

US State Codes

These are index IDs in the “State” input dataset (in the “Other” theme) -

This dataset contains standard FIPS codes and IDs for US States (see table below). These IDs are useful for many calculations where a single or list of US states has to be masked or processed.

Usage example: Population-2000 of New Hampshire where 33 is state ID code for NH -

(Pop00*CellArea)->where(State == 33)->sum

Table of US State FIPS codes and ID values-

River Basin codes

These are index IDs in the “RiverB” input dataset (in the “Other” theme) -

River (Watershed) Basins only with an area over 100,000 km2 are included in the table below (about 150 basins). Total number of basins in the RiverB dataset is over 10,000. So other basin IDs could be found by using i-Tool over the map above the Calculator -

Usage example: Population-2000 in Indus Basin where 26 is this basin ID code -

(Pop00*CellArea)->where(RiverB == 26)->sum

Table of River Basin ID Values and other Information

Continent codes

 These are index IDs in the “Cont” input dataset (in the “Other” theme) -

This dataset contains IDs of continents and the World ocean.

Usage example: Population-2000 of Africa where 1 is its ID code -

(Pop00*CellArea)->where(Cont == 1)->sum

Table of continent ID values

Questions and comments are welcome - alex.proussevitch@unh.edu

Written by Alex Prusevich
Edited by Richard Lammers, Alexander Shiklomanov, and Stanley Glidden