Map algebra on raster datasets represented by compact data structures

The increase in the size of data repositories has forced the design of new computing paradigms to be able to process large volumes of data in a reasonable amount of time. One of them is in-memory computing, which advocates storing all the data in main memory to avoid the disk I/O bottleneck. Compression is one of the key technologies for this approach. For raster data, a compact data structure, called k 2 -raster, have been recently been proposed. It compresses raster maps while still supporting fast retrieval of a given datum or a portion of the data directly from the compressed data. k 2 -raster’s original work introduced several queries in which it was superior to competitors. However, to be used as the basis of an in-memory system for raster data, it is mandatory to demonstrate its efficiency when performing more complex operations such as the map algebra operators. In this work, we present the


INTRODUCTION
In this work, we focus on spatial information represented with a raster model. Raster datasets represent geographic information by dividing space into a finite grid of cells and assigning a value to each cell. This includes images -such as remotely sensed imagery-, engineering, modeling, representations of parameters of the land such as pollution levels, atmospheric pressure, rain precipitations, land elevation, vegetation indices and so forth. The widespread use of different types of devices in many different scenarios causes the amount of information in raster format to increase rapidly. For example, only considering the family of satellites Sentinel, 20 TB of data are being generated each day, and the remote sensing repositories have reached the Peta Byte/Zetta Byte (PB/ZB) scale. 1 Using the sequential application of these operators, it is possible to perform complex data analysis on raster datasets. 32 The analysis procedures are categorized as those applied to an entire raster dataset (e.g., spatial autocorrelation) and those concentrated on specific locations: cell, neighborhood, or zones.
Depending on the neighborhood cells involved in computing the operator or function, there are four types of operators and functions [33][34][35] : (a) Local. The value in each cell of the output raster only depends on the values of the cells at the same location of the input datasets. Local operators are arithmetic, logical, and statistical functions (e.g., minimum, maximum, or mean). (b) Focal. In this case, the value of each output cell is computed taking into account the values of the cells at the same location of the input rasters and their neighboring cells in a specific range. Examples of focal functions include the computation of typical statistical functions (e.g., mean, mode, minimum, and maximum) corresponding to the neighborhood of a cell. Other processing such as convolution, kernel and moving windows also uses focal operations. (c) Zonal. In this type of operation, two input rasters are needed. The first one is used to divide the space into zones; then, the values of the output cells within a zone are computed processing the cells of the second input raster within that zone. (d) Global. This is the application of a bulk change to all cells. For example, adding a scalar value to all cells.
In the present study, we present the following operations. As a local operation, we show the point-wise arithmetic operation, see Figure 1A. In the case of global operations, we present two operations: thresholding obtains a binary raster where each cell of the output raster has a 1 if the value of the cell at the same position of the input raster is greater or equal than a given scalar value and a 0 otherwise (see Figure 1B), and scalar operation, where the output raster is the result of applying the same operation over all cells, such as adding a scalar value to all cells of the input raster (see Figure 1C). Finally, we present one zonal operation, zonal sum, which requires two input rasters, the zonal matrix and input matrix. The zonal matrix defines zones, which are formed by all cells having the same value. For each zone, the algorithm sums the cells of the input matrix, and then, the output is the zonal matrix, changing the original value by the sum, see Figure 1D.

Storage methods for rasters
A variety of formats and tools are available for the storage of raster data. Observe that a raster dataset can be easily stored as a matrix of values in generic matrix-processing tools, or even using raster-based image file formats. For example, GeoTIFF, a file format specifically designed for storing rasters, relies on the classical image file format TIFF enhanced with metadata to store geographical attributes of the raster. Next, we will introduce NetCDF, possibly the most widely-used file format for the representation of raster data.

NetCDF
Network Common Data Form (NetCDF) 36 is an Open Geospatial Consortium standard for creating, accessing, and sharing array-oriented scientific data, that is, rasters datasets. The data format is coupled with a set of software libraries that provide the ability to query the data. It is supported by most GIS software tools that handle raster data.
The file format arranges the cells of the raster in simple N-dimensional arrays. NetCDF provides a compressed or uncompressed representation of the data. Compression is based on Deflate, 37 which provides ten compression levels, ranging from level 0 (no compression) to level 9 (maximum compression). The main novelty of NetCDF is that even when compression is used, the data can be transparently accessed without performing an explicit decompression, that is, the access procedure is the same whether the data is compressed or not. Therefore, a trade-off arises between time and space depending on the used compression level; the more compression is chosen, the slower the querying process and vice versa. Moreover, queries over compressed data are efficient due to the use of a technique called chunking: data is compressed in blocks so that when a specific region of the data is demanded, only the relevant chunks need to be decompressed.

3.2
Compact data structures for raster data

Rank and select
Most of the compact data structures use two basic operations over bitmaps: rank b (B, p) counts the number of occurrences of bit b in bitmap B until position p and select b (B, n) returns the position of the nth occurrence of bit b in B. These operations can be efficiently solved in constant time 38 using n + o(n) bits of total space. In practice, only rank is solved in constant time whereas select is solved in O(log log n) time 39 with approximately 5% extra space over the original bitmap.

3.2.2
Representing rasters with k 2 -trees: k 2 -acc and k 3 -tree The k 2 -tree was initially designed to store and query web graphs. 15 However, actually, it is a region quadtree 40 for binary matrices built with the latest developments in the field of compact data structures.
Given a binary matrix M of size n × n, and a parameter k, M is represented with a tree. To build that tree, M is divided into k 2 submatrices of size n∕k × n∕k. A child node is added to the root node for each of those submatrices. Each node is labeled with one bit. A 0 indicates that the submatrix is full of 0-bits and a 1 means that the submatrix contains at least one 1-bit. Next, for each of the nodes labeled with a 1-bit, the process continues recursively.
To save space, the resulting tree is stored with two bitmaps without pointers, by following a level-wise traversal of the tree nodes: L is a bitmap formed by the bits corresponding to the last level of the tree, and T is a bitmap that contains the rest.
The k 2 -tree is limited to binary matrices. Thus, one way to represent rasters having values in the range v 1 < v 2 < · · · < v t is to use a k 2 -tree for each value. The representation is formed by t k 2 -trees K 1 , K 2 , … , K t , where each K i has a value 1 in the cells whose value is v ≤ v i in the original raster. This approach is the k 2 -acc. 21,22 Another way to represent a raster is to add a third dimension to the k 2 -tree. The k 3 -tree 21,22 stores points ⟨x, y, z⟩, where the first two values represent the position in the 2D space, and the third component is the value stored in that cell.

k 2 -raster
While the k 2 -acc and the k 3 -tree obtained better results both in space and search performance than other ways of representing rasters, such as those based on compressing GeoTIFF images, the k 2 -raster 23,25 has been proven superior to them, obtaining better results both in space and query times. The k 2 -raster follows the same quadtree-like building process of the k 2 -tree. Given a matrix M of integers of size n × n, † it builds a tree as follows. Given a parameter k, it divides M in k 2 submatrices of size n∕k × n∕k. Each submatrix adds a child to the root node of the tree storing the maximum and minimum values in the corresponding submatrix. If the maximum and the minimum values are equal, then the node becomes a leaf; otherwise, the process continues recursively on that submatrix.
Again, this tree is compactly represented using binary bitmaps. In addition, we also have sequences of integers corresponding to the maximum and minimum values, which will be also compressed.
We show in Figure 2 an example of k 2 -raster. In the upper part, we can see the matrix and its recursive subdivision. Just below, we can see the corresponding conceptual tree. The root represents the complete 8 × 8 matrix, and, as seen, the root contains the maximum and minimum values in that matrix. Under the label "Step 2," we can see the first division into four matrices of size 4 × 4, each producing a child of the root node. We highlight the maximum and minimum values in each submatrix, which are stored at the children of the root of the tree. The process continues until finding a submatrix filled with the same value or until we reach individual cells, which produce nodes with only the value of the cell (last level of the tree).
To reduce the magnitude of the numbers stored at the nodes, and thus open the possibility of representing them with fewer bits, the values at the nodes are stored as the difference with the same value at the parent. This is illustrated in the tree at center-bottom of Figure 2.
The topology of the tree is represented with a bitmap. Each node with children is represented with a 1-bit, whereas the leaves in levels before the last one are represented with a 0-bit. Thus, the bitmap is the simple level-wise traversal of the conceptual tree. Indeed, this is a simplified variant of LOUDS (level-ordered unary degree sequence) tree representation, 41 which is a compact representation for trees. In Figure 2, we can see this bitmap labeled as "Tree (T)." Using T, we can simulate the tree navigation with rank and select operations. More precisely, given a 1-bit at position p in T, that is, a node with children, these children are sequentially located from position children(p) = rank 1 (T, p) × k 2 of T, unless the children are individual cells, which are not represented in T. The parent of a node at position p of T is computed as parent(p) = select 1 (T, ⌊p∕k 2 ⌋). Due to the fast response time of rank and select, it is possible to efficiently access the value of a single cell, retrieve entire row/columns, or solve spatial range queries.
The values at the nodes are stored as arrays, following the same level-wise order of the tree. We need two arrays Lmax for maximum values and Lmin for the minimum. These arrays are compressed with Directly Addressable Codes (DACs), 19 a compression scheme for integer sequences that provides the ability to direct access to any given position. We use the † If the matrix is not a square, we can fill it with "empty cells" until it becomes a square. This does not impact the space requirements of the final representation, as the k 2 -raster is very efficient representing regions with the same value.
DACs version that calculates the optimal number of bits at each level. DACs exploit the fact that encoded differences tend to be small. This is shown at the bottom of Figure 2, as part of the final representation.
As explained, the k 2 -raster obtains a very compressed representation of the raster matrix and, at the same time, it indexes the data, enabling fast queries over the raster matrix. It has proven to be superior to the rest of the state-of-the-art techniques both in compression power and query performance, except in the case of NetCDF's compression power, since they are almost on par. 23 Compared to k 2 -acc, the technique used in the previous proposal of map algebra using compact data structures, 31 the k 2 -raster not only obtains less space consumption and query times, but it also scales better when increasing the size of the input data or the number of different values. In fact, when the number of different values is moderate, the current implementation of k 2 -acc does not work, and any possible implementation will have problems to deal with that situation since it requires a complete k 2 -tree for each different cell value present in the raster. Moreover, except in rasters with a F I G U R E 2 k 2 -raster example. Example of integer raster matrix (top), conceptual tree of the k 2 -raster (center-top), conceptual tree using differential encoding (center-bottom), and final representation of the raster matrix using compact data structures (bottom). rMax and rMin denote the maximum and minimum values of the root node. Lmax and Lmin contain the maximum and minimum values of each node, following a level-wise order and using differential encoding. This example uses k = 2. very low number of different values in the raster, k 2 -acc obtains between 2.5 and 8.75 times worse compression. In query times, in most queries, k 2 -raster is faster than k 2 -acc, up to one order of magnitude.

Quadtrees
k 2 -acc, k 3 -tree, and k 2 -raster use a quadtree as the basic data structure (octree in the case of k 3 -tree), which is then enhanced with several other components. Quadtrees were originally introduced to compress images. 42,43 However, since then, many different variants have appeared. We can roughly classify them in those focused on indexing data [44][45][46] and those focused on compressing images or 3D meshes. [47][48][49][50][51][52][53] As indexes, typically, data are stored in an array-like data structure and a quadtree is used as an auxiliary index, 44,54 whereas, in the case of compression, no indexation is applied. However, k 2 -acc, k 3 -tree, and k 2 -raster combine in the same data structure two indexes and the storage of the data. A similar fusion has been proposed more recently, 55 however this system only indexes the space, whereas k 2 -acc, k 3 -tree, and k 2 -raster, in addition to the spatial index (quadtree), are equipped with an index on the values at cells.
Focusing on compression, the basic idea of the quadtree is to take advantage of the intrinsic data smoothness, by representing parts of the image with only one value or, at least, represent parts of image with less information than the original by taking advance of the redundancy between close parts of the image. k 2 -raster uses the same idea, but adding some other compression mechanisms.

3.4
Operations over rasters using compact data structures Brisaboa et al. 56 proposed algorithms to efficiently perform Boolean operations over binary raster datasets represented using k 2 -trees. They included algorithms for computing the union, intersection, difference, and complement of binary rasters using two different variants of the k 2 -tree representation for binary relations. Thus, it is related to our work as they propose algorithms for computing some local operations over raster datasets, but it differs from us in the sense that their work is focused on binary rasters, whereas we extend the applicability of these types of map algebra operators for rasters of any nature. Later, Quijada-Fuentes et al. 57 presented algorithms for solving set operations over k 2 -tree and over another compact data structure called binary relation wavelet tree (BRWT); however, again, this work only considers binary rasters. The closest work to ours is the development of map algebra operations on rasters represented using k 2 -acc. 31 They presented algorithms for the following operations: thresholding, sum and multiplication by a scalar, point-wise sum, and zonal sum. However, as explained, k 2 -acc has poor performance in comparison with k 2 -raster, especially when the raster has a moderate or high number of different values, which is the most common situation.

OUR PROPOSAL
First of all, let us introduce the following notation. Let p M 1 be a pointer to a node of the k 2 -raster M 1 . With * (p M 1 ), we access the pointed node. That is, * (p M 1 ).min and * (p M 1 ).max return the minimum and maximum values, respectively, stored at the pointed node.

Global operations: Arithmetic operation by a scalar value
This operation performs the addition, subtraction, division, or multiplication of a scalar value to all cells of the matrix. It is performed differently depending on whether the operation is addition/subtraction or division/multiplication. However, in both cases, it is simple.
In the case of addition or subtraction, we only have to sum or subtract the input scalar to the values maximum and minimum at the root node, that is, the fields rMax and rMin of the Figure 2. Observe that, independently of the size of the raster, this implies just two operations.
In the case of division and multiplication, the algorithm requires applying the operation to all nodes of the tree. This implies applying the operation to rMax and rMin and to all the values in Lmax and Lmin. However, the structure of the tree does not change; thus, the operations over Lmax and Lmin can be done sequentially without traversing the tree.

Local operations: Point-wise
In this case, the value of each cell of the output matrix M o (i, j) is computed as the operation over the cells at the same location of the two input matrices, that is, , being Θ an arithmetic operation (+, −, *, or /). Algorithm 1 shows the pseudocode of the operation. At any given step of the algorithm, it manages two pointers. Initially, both pointers point to the root node of the two input k 2 -rasters. The algorithm traverses downwards both input trees until reaching quadrants of size 1 × 1, that is, individual cells, or a uniform quadrant. It is invoked as Point-wise(Θ, n, 1, p M 1 , p M 2 ). The first parameter is the operation, n is the (possibly extended) raster matrix size, the third parameter is the current level, and p M 1 and p M 2 are pointers to the current node of each input k 2 -raster. In the initial call, both point to the root node of the corresponding k 2 -raster.
k, T , Vmax , and Vmin are global variables. k is the k value of the trees. T , Vmax , and Vmin are lists, initially empty; there are T , Vmax , and Vmin lists for each level of the output tree. In addition, the global variables pmax and pmin indicate the last written position of Vmax and Vmin , respectively.
After running the algorithm, all T sequences must be joined to make up T. The same must be done with Vmax and Vmin to obtain Vmax and Vmin, which, in turn, must be converted into Lmax and Lmin by computing the differences and then encoded using DACs. Observe that the algorithm returns the maximum and minimum values of the resulting matrix, that is, rMax and rMin.
Each call to Point-wise processes one quadrant (let us say q p ) of both input matrices, initially the whole raster, to obtain the same quadrant of the result.
If the two input pointers correspond to uniform quadrants, that is, all cells have the same value, then lines 4 and 5 obtain the value that will have the output raster in all cells of the processed quadrant. Line 6 performs the operation and introduces the resulting value in both maxval and minval (the answer of the call). When at least one of the nodes is not a leaf, for each non-leaf node, the fors of lines 8 and 9 process its k 2 children, with lines 10-17 obtaining the pointers to those children. In the case of leaf node, the pointer remains pointing to the current processed node q p . For each subquadrant of the processed node, line 18 performs a recursive call to process it.
After the recursive call, line 19 adds the returned maximum value to Vmax . Line 20 checks if the returned values are different. In that case, it appends the minimum value to Vmin and sets up a 1 in the T of that level. If the maximum and minimum values are equal, it sets up a 0 in T , provided that we are not returning from a recursive call that processed an individual cell. Lines 28-30 keep track of the minimum and maximum values found so far during the computation of q p of the target raster.
After processing the k 2 children (q p1 , q p2 , … , q pk 2 ), line 34 checks whether the maximum and minimum values of all children are equal, which indicates that all the children contain the same value. Thus, the algorithm must undo the last operations, as these nodes will not have a representation in the data structure. This can be easily done by removing the last k 2 positions of T and Vmax , or just moving the pointer that indicates their last written position, k 2 positions backwards (line 35). Finally, the algorithm returns the maximum and minimum values to its parent call. Figure 3 shows two rasters, the output of their point-wise sum and their conceptual k 2 -rasters. We are going to use them to illustrate how Algorithm 1 works. The first level division of the rasters is depicted with thick lines. The second level division is drawn with thin lines. First level division produces four quadrants labeled as q 1 , q 2 , q 3 , and q 4 .
In order to distinguish between the two input rasters, we will use q i (M 1 ) to refer to the quadrant q i of the first input raster and q i (M 2 ) for the other. To refer the corresponding node of quadrant q i (M j ), we will use n i (M j ). Equally, where a pointer is required, by n i (M j ) also stands as the pointer to that node in the k 2 -raster.
At the returning of the call, Vmax 1 [1] receives the 5 and line 25 sets T 1 [1] to 0 indicating that it is a leaf node (see step 3 of Table 1). Lines 28 and 30 set minval and maxval to 5, since that is the minimum and maximum values found so far in the children of the root.
We returned to the call of step 1, and now the for of line 10 sets j to 1, and this implies that the subquadrant 2 must be processed, and thus, the recursive call Point-wise(+, 2, 2, n 2 (M 1 ), n 2 (M 2 )) is issued (see step 4 of Table 1). Since q 2 is not uniform in both input trees, the flow reaches line 8, and then the k 2 subquadrants of q 2 must be processed with recursive calls. Steps 5-12 of Table 1 show the k 2 recursive calls for processing the subquadrants. Step Call Actions Rec call Actions Rec call 12 Vmax 2 (9,12,11,8) 13 18 Step 5 shows the call Point-wise(+, 1, 3, n 11 (M 1 ), n 11 (M 2 )). n 11 is an individual cell, with value 2 in M 1 and 7 in M 2 , therefore lines 4-6 set maxval and minval to 9. After the recursive call, when the flow returns to the call of step 4, line 19 adds the 9 to Vmax 2 (see step 6). Neither Vmin 2 nor T 2 are updated since we are at the last level (individual cells). Until step 12, the same procedure is applied to the rest of the children of n 2 .
After processing all children of n 2 , the flow returns to the call of step 1 in line 19. The call Point-wise(+, 2, 2, n 2 (M 1 ), n 2 (M 2 )) returns (12,8), and thus 12 is attached to Vmax 1 and 8 to Vmin 1 . Finally, a 1 is attached to T 1 since n 2 in the output tree has children.
Finally, the process of q 4 is similar to the process of q 2 since, although in M 1 that node is a leaf, it is not in the case of M 2 (see steps 16-25).

Global operations: Thresholding
Algorithm 2 shows the thresholding operation. It is also a recursive procedure that starts at the root but, in this case, it only traverses one tree. It is invoked as Thresholdin(n, , p M , t). The first parameter is n the raster matrix size, the second parameter is the current level, p M is a pointer to the current node of the input k 2 -raster, and t is the threshold. In the initial call, it points to the root node. The global variables are the same as in the case of the Algorithm 1.
The main difference is the recursion stop condition. As in the case of the point-wise operation, when we reach a leaf, recursion stops, but now returning 0 or 1 in both maxval and minval depending on whether the value at the node is greater or lower than the threshold.
There are two additional cases where recursion stops. If the maximum value at the processed node is smaller than the threshold, then it is sure that all values of that subquadrant are smaller than the threshold, so the recursion stops returning a 0 in both maxval and minval. Observe that, since we store the same value in both maxval and minval, then at the returning point of the recursion lines 18-35 will create a leaf node. The other case is when the minimum value of the processed node is greater than the threshold, which results in returning a 1 in both maxval and minval. Otherwise, recursion continues performing a call for each of the k 2 children of the processed node (lines [14][15][16][17], and after returning from that call, the output tree is built (lines [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34].
Note that this operation takes advantage of the k 2 -raster index capability provided by the min-max values at nodes, which avoids processing tree branches that can be solved at the upper levels of the tree.

Zonal operations: Zonal sum
Observe that the point-wise and thresholding operations use recursive algorithms that start at the root of the input trees, and process them downwards until they reach the leaves. Then, when returning from the recursive calls, they build the output tree from the leaves until the root. Therefore, basically, the process involves a downward and an upward traversal.
In this operation, this approach cannot be used, because we have to sum all values in the input matrix overlapping any given value of the zonal matrix. Therefore, until processing completely both input matrices, we cannot start the construction of the output tree.
Zonal sum will be solved with two traversals as well, the first one traverses the input trees downwards and the second one builds the output tree upwards. However, the first traversal must finish completely before the second one starts. This prevents the use of a recursive procedure that mixes the two traversals. The first traversal to obtain the zonal sum of each value of the zonal matrix uses a recursive algorithm. The second traversal is a sequential process of the levels of the zonal matrix from the leaves until the root to build the output tree.
Note that the output tree is basically the zonal tree, except that some non-uniform submatrices (non-leaf nodes) of the zonal tree may become uniform, and thus leaves, in the output tree. Therefore, during the first downward traversal, in addition to obtaining the zones, the algorithm records the structure of the zonal tree using node lists. Each list contains the nodes of one level from left to right. Then, the upward traversal of the zonal matrix is performed using those lists. maxval, minval ← 1 /*All values in the quadrant are above t*/ 13: else 14: for i ← 0 … k − 1 do 15: for j ← 0 … k − 1 do 16: p M ← Child(p M , i ⋅ k + j) /*Obtains the child*/ As explained, the upward traversal replicates the zonal tree, and only when the algorithm finds that a non-uniform quadrant of the zonal matrix becomes uniform in the output, then it changes the output accordingly. To detect those situations, the lists, in addition to informing about the structure of the zonal tree, are also used as an auxiliary structure during the computation of the output to store the min-max values of the nodes of the output tree.
Algorithm 3 shows the pseudocode. The hash table Hash has all its entries initialized to 0. Vmax2 is the Vmax array of the zonal matrix M 2 , and it is a global variable. The lists L 1 , L 2 , … , and P 1 , P 2 , … are also global variables. As in the previous algorithm, T , Vmax , and Vmin are lists, initially empty; there are T , Vmax , and Vmin lists for each level of the output tree. Again, after running the algorithm, all T sequences must be joined to make up T and Vmax and Vmin must be processed to obtain Lmax and Lmin. end if 25: if minval>min thenminval ← min 26: end if 27: if (i + 1) mod k 2 = 0 then /*The parent must be updated*/
The algorithm is called as ZonalSum(n, # , p M 1 , p M 2 ). The first parameter is the raster matrix size, # is the number of levels of the zonal tree (M 2 ). p M 1 and p M 2 are pointers to the root node of both input k 2 -rasters.
Line 1 performs the call to the recursive process that adds the values of the input matrix for each different value in the zonal matrix. The returned hash table Hash contains the sum for each value of the zonal matrix.
Algorithm 4 shows the pseudocode of this process. It is basically the same recursive process of the point-wise operation, excluding the construction of the answer. The differences are two. First, in line 4, when the recursive process reaches a leaf in both trees, the hash entry of the value of the zonal matrix accumulates the values of the input matrix. Second, when a node of the zonal matrix is accessed, line 17 adds it as an entry of the list of nodes L +1 . The entry contains a pointer to that node and minimum and maximum values set to -1, and lines 14 and 15 add the position of the node within its level (from left to right), if it is a parent node. This list of parent positions is required during the construction process of the output.
After that call, returning to Algorithm 3, the while of line 4 processes the levels of the zonal matrix from the last level until level 2 (that of the children of the root). Variables pmax and pmin are the positions of Vmax and Vmin where a new element must be written. The while of line 8 processes sequentially the entries of L .
If the processed entry of L has a -1 in the max field, then that means that it is a leaf in the zonal matrix; otherwise, its children would have been processed during the level + 1 traversal, and this would update the max and min fields of the processed entry. Therefore, in line 10, the corresponding position in Vmax is filled with the sum associated with the value of the zonal matrix. Otherwise, the maximum and minimum values are already at the entry, and thus, Vmax and Vmin are updated accordingly (lines 12-16).
Next, in lines 18-22, T is updated. Lines 23-26 keep track of the maximum and minimum values of the children of the current parent node. Line 27 checks when the traversal reaches the last child of a parent node. In that case, we have to update the minimum and maximum values at the entry of the parent node in the upper level. Before that, if we find that the maximum and minimum values of the children are the same, line 29 moves the pointer pmax k 2 positions backwards to erase the children of the output tree, since the parent node will be a leaf.
Lines 31 and 32 obtain the position in level − 1 of the parent of the processed children using the list where we stored the position in each level of the parent nodes. Lines 33 and 34 update the maximum and minimum values at the parent node.
Next, we are going to illustrate the Algorithm 3 using the rasters of Figure 4. The algorithm is invoked as ZonalSum(4, 3, p root (M 1 ), p root (M 2 )). Line 1 issues Sum(4, 1, p root (M 1 ), p root (M 2 )). As seen, this is a recursive process that obtains the hash table with the sum for each value of the zonal matrix, the lists L having the structure of the zonal tree by levels, and the lists P of positions in each level having parent nodes. Figure 5 shows these elements for our example. Observe that the lists L store the positions in Vmax2 corresponding to the nodes of the level and that the min and max values are set to -1. The list of parents is our example is only available for level 2, and, as seen, it marks only positions 1 and 3 (starting at position 0), since in the k 2 -tree of the zonal matrix, those positions of level 2 contain parent nodes (the nodes labeled 9-6 and 18-10).
After that call, returning to Algorithm 3, the while of line 4 processes the levels of the zonal tree starting at the last one (3). The while of line 8 processes L 3 from left to right. The first entry is ⟨5, −1, −1⟩; since it has a -1 in max, the algorithm obtains from the hash table the sum associated with the zonal value. For this, first, the zonal value is obtained using the pointer, so the zonal value is Vmax2 [5] = 7. Next, the sum associated with 7 is obtained from the hash table, which in our case is a 2 (see the second entry from the top of table Hash of Figure 5). Finally, that sum is added to Vmax 3 . Lines 18-22 are not executed since we are at the last level and it is not represented in T (see row with i = 0 and = 3 of Table 2). Next, maxval and minval, which keep track of the maximum and minimum values of the parent node of the currently processed node, are set to 2. The next entry is ⟨6, −1, −1⟩, which adds a 3 to Vmax 3 , and updates maxval to 3. The same process is carried out with the following two entries of L 3 , then when reaching i = 3, this results in Vmax 3 = (2, 3, 3, 2) (see rows with i = 1, 2, 3 and = 3 of Table 2). However, at that point, that is, when processing the entry ⟨8, −1, −1⟩, we are processing the last child of the parent node labeled 9-6 in the zonal tree. Therefore  Table 2. See in Table 2 that from i = 4 to i = 7 for = 3, the k 2 leaves of the node 18-10 of the zonal tree (which are the values 11, 12, 18, 10) have a zonal sum of 1, and thus the four 1s are added to Vmax 3 . When processing the fourth child, that is when i = 7, then the if of line 27 becomes true; minval = maxval = 1, and thus, pmax is moved k 2 positions backwards, thus removing the last four added 1s. In addition, lines 31-34 update the minimum and maximum values of the last entry of L 2 , setting both of them to 1.
Next, the flow returns to line 4 to process level 2. The first entry of L 2 (⟨1, −1, −1⟩) has a -1 in max, so it is a leaf. Using the pointer, we obtain the associated sum (12) from the hash table (see entry with = 2 and i = 0 of Table 2). However, now we are not at the last level; thus, lines 18-22 are executed, adding in this case a 0 to T 2 , which indicates that this node is a leaf. The next entry is ⟨2, 3, 2⟩. Given that max is not -1, the node is not a leaf, and then Vmax 2 and Vmin 2 are updated with the values stored at the entry, as it can be seen in the row with i = 1 and = 2. In addition, a 1 is added to T 2 indicating that it is not a leaf. The third entry is like the first one. The fourth entry (⟨4, 1, 1⟩) has max different from -1. In this case, both values are the same (1) and thus only Vmax 2 and T 2 are updated (see row with i = 3 and = 2 of Table 2).

A note on complexity
In this section, we will analyze the space and time complexity of the proposed algorithms. In the case of the space complexity, we do not include the input/output sizes, but the additional space required to compute the map algebra operations and build the output.

Global operations: Arithmetic operation by a scalar value
In the case of adding or subtracting a scalar value, both space and time complexities are O(1), as only the fields rMax and rMin need to be modified.
In the case of multiplication by a scalar value, not only rMax and rMin are modified, by all the values of L max and L min arrays. Thus, time complexity is linear to the time required to build these arrays, that is, O(|L max | + |L min |), where L max and L min are the arrays storing the maximum and minimum values of the input k 2 -raster. Regarding space, it requires an additional space for building the DACs arrays of the output, as they cannot be created in compact space. Thus, the space complexity is O(|L out max | + |L out min |), being L out max and L out min the arrays storing the maximum and minimum values of the output k 2 -raster.

4.5.2
Local operations: Point-wise The time complexity of the algorithm is linear to the input size, as it just process once each node of the input k 2 -rasters. Thus, it is O(|M 1 | + |M 2 |), being M 1 and M 2 the input k 2 -rasters. Regarding space complexity, it just requires the additional space for building the DACs arrays of the output, that is, O(|L out max | + |L out min |), being L out max and L out min the arrays storing the maximum and minimum values of the output k 2 -raster.

Global operations: Thresholding
The time complexity of the algorithm is linear to the input size, as it just process once each node of the input. Thus, it is O(|M|), being M the input k 2 -raster. Regarding space complexity, it just requires an additional space O(|L out max | + |L out min |) for building the DACs arrays of the output, being L out max and L out min the arrays storing the maximum and minimum values of the output k 2 -raster.

Zonal operations: Zonal sum
Again, the time complexity of the algorithm is linear to the input size, as it just process once each node of the input. Thus, it is O(|M 1 | + |M 2 |), being M 1 the input k 2 -raster and M 2 the zonal matrix. Assuming M 2 is smaller than M 1 , time complexity is O(|M 1 |).
Regarding space complexity, again it just requires an additional space O(|L out max | + |L out min |) for building the DACs arrays of the output, being L out max and L out min the arrays storing the maximum and minimum values of the output k 2 -raster.

EXPERIMENTAL EVALUATION
We include in this section to three different experiments. In the first one, we compare our algorithms with two baselines: NetCDF and a naive solution that completely decompresses the input rasters, runs the query over the decompressed data, and finally compresses the result to obtain a k 2 -raster. The second experiment compares our method with the algorithms proposed for the k 2 -acc. 31 The reason of this division is that we use the authors' implementation of the algorithms on k 2 -acc, which do not write the output, that is, they only compute the result in main memory, but they do not write that result back to disk. Since NetCDF does that, in order to be fair, we decided to run a second experiment with our algorithms modified to behave in that way. Finally, the third experiment compares the memory consumption of all methods jointly, since the process of writing to disk does not affect this parameter.
In the case of netCDF, we used the netCDF library ‡ (v.4.7.4) with the recommended level 2 of compression, which obtains a good trade-off between space and access times. For the second experiment, the source code for the algorithms of the map algebra operations over the k 2 -acc was provided by their authors. 31 It uses k = 2.
Our algorithms were implemented in C++ using the SDSL § library. 58 All the experiments were run on a dedicated Intel® Core TM i7-3820 CPU @ 3.60GHz (4 cores) with cache sizes 32KB (L1), 256KB (L2), and 10MB (L3), and 64GB of RAM. The operating system was Debian 9.12 with kernel 4.9.0-9-amd64. The code was compiled with gcc version 6.3.0 with -O3. We measured elapsed time (in seconds) and, in our experiments, all data were initially on disk. For all experiments, we repeat the map operation 30 times and report average times. Input datasets and values for scalar and thresholding operators were chosen randomly.

Datasets
In our experiments, we used real datasets from two sources: • WorldClim ¶ dataset 59 provides several layers of global climate information. The surface of the world is divided into equal-sized tiles. In turn, each tile is a raster with cells covering about 1 square kilometer. In our case, we used the layer with mean temperature. Whereas in the original dataset it is represented by a real number with one decimal, in our experiments we transformed it into an integer by multiplying the value by 10. The collections obtained from this source are denoted clim.
• Spanish Geographic Institute # (SGI) provides a DTM (Digital Terrain Model) of Spain, that is, spatial elevation data of the terrain. The surface is divided into rectangular equal-spaced tiles with 5 m of spatial resolution. Each cell of a tile stores a real number of at most three decimal digits.
Since the presence of uniform areas, that is, areas with the same value, may benefit the k 2 -raster, we considered two different levels of precision for the datasets from the second source by using different numbers of decimal digits. We generate collection dfm 0 by truncating the decimal digits of the raster values, and collection dfm 3 by using the original numbers, but represented as integers, that is, using three decimal digits and multiplying the original value by 1000.
We have created several collections of rasters. In each collection, all rasters have the same size. We decided to use collections of several rasters in order to eliminate any bias due to the use of a single matrix. Our collections are of different sizes to analyze the scalability of the algorithms. In collection labeled as 1x1, each raster is a matrix built using just 1 tile of the original dataset, 2x2 collections contain matrices built using 2×2 adjacent tiles, and so on. Tables 3 and 4 show the characteristics of our collections of rasters. The data shown in the tables represent the mean values of all rasters in each collection. We also include in these tables the space consumption (in MBs) required by each of the three methods used, that is, k 2 -raster, NetCDF and k 2 -acc when representing these collections. We do not include the results of k 2 -acc for dfm 3 collection since it has too many different values and this technique is not able to run over these datasets. Figure 6A shows the times of the sum by a scalar. This operation benefits clearly to k 2 -raster since the algorithm only implies adding (or subtracting) the scalar to the maximum and minimum values of the root node, which is a constant time algorithm. However, observe that the slope of the line corresponding to k 2 -raster slightly increases as the size of the dataset also increases. This is because our algorithm loads the complete k 2 -raster into main memory, and thus as the size of the k 2 -raster increases, the time for loading and writing from/to disk also increases. Observe as the uniformity of the raster benefits in any case to the k 2 -raster. The most uniform collection is clim; thus, our algorithm is between 2.9 and 17.6 times faster than the naive approach and between 2.9 and 6.9 times faster than NetCDF. Something similar happens with dfm 0 , whereas with the least uniform, dfm 3 , improvements range between 4.1 and 6 times in the case of the naive approach and between 2.5 and 2.7 times in the case of NetCDF. Again, in this operation, the changes are due to the effect of uniformity in the compression power of k 2 -raster, as the more uniformity is found in the original raster, the shorter is the tree and thus more compression is achieved, and, as explained, in this operation, the main cost for k 2 -raster is to load and write from/to disk. Figure 6B shows the results of product by a scalar. In this case, all values in the leaves of the tree must be accessed and operated. Although the tree is not accessed, and only arrays Lmax and Lmin are sequentially accessed, in this case, k 2 -raster is not able to outperform NetCDF in most cases. On the other hand, it achieves better performance and scales better compared with the naive approach. Our algorithm is between 1.38 and 2.3 times faster than the naive approach, whereas compared with NetCDF, except in the smallest collection of clim where they are on a par, in the rest, it is between 1.05 and 1.38 times slower.

NetCDF experiment
In general, when a sequential scan of the whole raster is needed, NetCDF is capable of matching the speed of the k 2 -raster and even surpassing it. Several factors explain this: • k 2 -raster is an index also carrying the data. However, when all data must be accessed, the advantages provided by the indexes disappear.
• Once the indexes are useless, we have to pay attention to the speed of the underlying compression method. NetCDF uses Deflate, which is a fast decompressor, whereas the DACs of the k 2 -raster are slower. This is expected, as Deflate does not worry about direct access, since it is an archival method, and NetCDF has to compress by blocks in order to obtain some sort of direct access. However, DACs are precisely Directly Addressable Codes and any individual value can be decompressed separately.
• The third element that affects speed is that the k 2 -raster in an in-memory method, which favors k 2 -raster compared to NetCDF, a disk-based method. Figure 6C shows the results of the point-wise sum operation. This is precisely the best example of the previous explanation. Here, the indexes of the k 2 -raster are not helpful since the two input rasters must be completely accessed. However, the in-memory approach of k 2 -raster balances the comparison, and, at the end, there is no clear winner.
Thresholding is precisely the other side of the coin, as shown in Figure 6D. Here, the indexes are capable of processing big parts of the input raster in the upper parts of the tree with just one step. That is, if we know that all cells in a subtree are below the given threshold, the algorithm simply sets them all to zero, without accessing them. Therefore, our algorithm is between 2.1 and 17.6 times faster than the naive approach and between 1.63 and 6 times faster than NetCDF, since both have to sequentially process cell by cell the complete input raster.
In zonal sum, the uniformity has even more impact. Recall that, by construction, as the raster is more uniform, the corresponding k 2 -raster is smaller. Thus, since operations require a top-down traversal, the time performance improves. However, in this operation, the uniformity of the zonal raster is even more critical since, in addition to the aforementioned effect, the algorithm produces an output close to that raster, and the second part of the algorithm takes time proportional to the number of zones.
This can be seen in Figure 7A,B. Figure 7A has a zonal raster with only 10 zones; thus, in the most uniform collections, dfm 0 and clim, our algorithm is able to outperform NetCDF. More concretely, it is between 1.11 and 1.34 times faster in dfm 0 and between 1.87 and 2 times faster in clim. However, in dfm 3 , our algorithm is between 1.07 and 1.37 times slower. We can see the impact of increasing the number of zones to 500 in Figure 7B. Now, our algorithm is slower in dfm collections, between 1.31 and 1.51 times slower, and it is on a par in clim. The effect of the number of zones of the zonal raster can be seen in Figure 7C, where we show the time versus number of zones with the collections of size 4x4. In all cases, our algorithm outperforms the naive approach. With respect to NetCDF, again the uniformity is the key factor, and depending on it, one is better than the other.

5.3
k 2 -acc experiment k 2 -acc implementation has serious limitations to deal with datasets of moderate size, not only in the number of cells, but also in the number of different values in the dataset. Therefore, we only present results with the dfm 0 and clim datasets, as with dfm 3 , the k 2 -acc is not able to run or requires extremely long execution times. Figure 8A,B show the scalar operations. In these operations the k 2 -acc can compete with k 2 -raster. The k 2 -acc is based on having one k 2 -tree per value present in the raster, which marks the cells where that value is present. Therefore, for each different value, that value and a k 2 -tree are stored. To solve these operations, only the values are modified, while the k 2 -trees remain unchanged. Namely, the arithmetic operation is applied to the values, and thus the cost of the operation is limited to an arithmetic operation for each different value.
However, in the case of the sum, k 2 -raster only has to perform two sums, one for the minimum value and another one for the maximum. This is constant independently of the number of different values in the raster; thus, as seen, k 2 -raster is around two-three times faster, except in the smallest dataset.
Nevertheless, in the case of the product, k 2 -raster has to traverse all the leaves, and thus, k 2 -acc is between 1.74 and 8.71 times faster.
However, as soon as the k 2 -trees must be accessed, k 2 -acc becomes extremely slow. In the case of the point-wise operations, Table 5 shows the results only for the datasets of size 1 × 1, as with bigger datasets, k 2 -acc did not run. As seen, differences are of at least two orders of magnitude.
Thresholding is the other operation where the k 2 -acc can compete. In fact, it is very easy to solve it, as it only requires to recover the k 2 -tree representing the input threshold. Still, k 2 -raster is faster, as shown in Figure 8C, with the exception of the smallest datasets. Observe that k 2 -tree traverses the input tree with the help of the minimum and maximum values in the nodes, which allow the fill of large zones of the output tree with 0s or 1s in the upper levels of the tree; thus, the output tree will probably be a small tree. Still, this implies more processing that a simple copy of the k 2 -tree corresponding to the threshold value. However, the k 2 -acc is harmed by its poor compression ratios, as seen in Tables 3 and 4, which implies slower load times in main memory than in the case of k 2 -raster.
In the case of the zonal sum, we only include experiments with 10 and 100 zones, as the k 2 -acc requires large running times. Again, since the values of the cells must be accessed, k 2 -acc has serious problems. Figure 9A,B show values where k 2 -raster can reach improvements of up to two orders of magnitude.

Memory consumption
In this experiment, we considered resident memory. For memory consumption, the disk-based classical approach of loading one block at a time is the best scenario for NetCDF. In fact, the in-memory approaches try to avoid this traditional method derived from the low memory capacity of old computers. The decreasing in price of memories opens the opportunity of keeping all data in main memory all the time, thus taking advantage of the faster memory times, around 6 orders of magnitude faster. Nevertheless, compression is needed to be able to fit all data in main memory, which lowers that gap. The k 2 -raster's original paper shows astonishing improvements in query times when its indexes are useful, whereas in this work, it shows variable results when all data must be sequentially processed and the indexes are useless. However, as a consequence of the change of approach, regarding memory consumption, see in Figures 10 and 11 that NetCDF is unbeatable. Still, we show that the algorithms presented here are noticeable better that the naive approach, therefore if one chooses the k 2 -raster as its storage method, our new algorithms obtain big memory savings when comparing to the naive ones.
With respect to to the k 2 -acc, we do not include data from the point wise operation and, in the rest of operations, we do not include data of the dfm 3 dataset due to the problems of k 2 -acc running that operation and datasets. Unlike time, where the k 2 -raster is a landslide winner, in this parameter there is no clear winner.

CONCLUSIONS AND FUTURE WORK
With compression power similar to that of NetCDF, k 2 -raster has been shown clearly superior when issuing the most common queries, namely, retrieving a region of the space and retrieving the cells in a region of the space with values in a given range. However, although less common, the capability for efficiently solving map algebra operations must be investigated to back the hypothesis of using the k 2 -raster as the method for representing rasters all the time, and more precisely in in-memory environments.
In this work, we have proposed efficient algorithms for performing map algebra operations over the compressed representation of the input rasters using k 2 -raster, and proved, in the experimental evaluation, that our algorithms are preferable to the naive approach of completely decompressing the rasters, running the operation with the uncompressed rasters, and then compressing the result, which is the key of an in-memory system. In addition, to complete the experimental evaluation, we have also included the classical method to store rasters, that is, NetCDF and a previous compact data structure for storing rasters, k 2 -acc.
In the case of NetCDF, there is no overall winner, as there are operations where k 2 -raster obtains the best performance, and there are other operations for which NetCDF outperforms the proposed algorithms using k 2 -rasters.
The proposed algorithms using k 2 -rasters obtain the best performance for summing/subtracting by a scalar and thresholding. In the first case, the k 2 -raster benefits from its structural properties and how it compactly represents the data, solving that operation in constant time. In the second case, the k 2 -raster takes advantages of its indexes. However, the proposed algorithms are not capable of clearly overcome NetCDF when the rasters must be processed sequentially, such as for the point-wise operation of zonal sum. In those cases, the algorithms do not use the indexes of the k 2 -rasters; thus, traversing their tree shape is less efficient than processing the array-like structure of the NetCDF, and only the in-memory approach of k 2 -raster is able to balance the times.
In the case of k 2 -acc, k 2 -raster is the clear winner in time. With much better compression power, it only loses in one operation, the product by a scalar, but in those operations where the k 2 -acc has to determine the value of the cells, k 2 -raster is orders of magnitude faster.
We did not include focal operations, as these operations can be solved with the window query shown in the k 2 -raster's original paper. When the focal operation only affects to one cell, such as computing the slope of the terrain around a cell, that operation will be very efficient, since the k 2 -raster is very fast extracting a portion of the data (window query) given that it is in essence a quadtree, that is, a spatial index. However, when the operation affects to all cells of the raster, such as computing the maximum, minimum, average and so forth of all cells within an area around each cell, the best alternative is to apply the window query to the complete raster, as it does not make any sense to issue a different window query for each cell. This is precisely the naive baseline presented in our experiments, which can be used in any operation, as those cases where a native k 2 -raster algorithm would be inefficient.
As a summary, putting in a balance the pros and the cons, k 2 -raster remains as a better choice: it obtains almost the same compression performance as NetCDF and it is much better when issuing common queries, 23,25 other complex operations, 26 and, as proved in this work, also some of the map algebra operations, and when it is not the clear winner, the performance is similar to that of NetCDF. As future work, we plan to analyze the parallelization of the proposed map algebra implementation, as many of these algorithms are clearly parallelizable.