The pros and cons of the GMT grid padding

The GMT team will be discussing the perceived problems of the GMT padded grids, especially in the context of external environments passing grids in and out of GMT. This issue tries to explain what the pad is, why it is there, what it buys us, and what the concerns may be, as well as solutions.

**What is the grid padding?**

Before the birth of the external tools built on top of GMT (GMT/MEX, PyGMT, Julia), we made the decision that internally, GMT grids should be stored in memory as a 2-D matrix that is extended by two extra rows at the top and bottom, and two extra columns at the left and right sides. Thus, a data grid of size _ny_ by _nx_ is stored as a matrix of size _my_ by _mx_, where typically _mx = nx + 4_ and _my = ny + 4_. Generally designed to deal with files, GMT will read in a file (which of course has no padding) and place it into the padded array. Once it is time for writing the output to a new file, we only write out the original size and the pad is never written to file.

**Why is there a pad?**

In one word, coding _simplicity_. Many GMT grid operations, but not all, involve algorithms that require a grid to have boundary conditions of second order. This means the grid _z(x,y)_ has to meet two conditions on each side. These boundary conditions can be translated into the values of the two artificial nodes beyond the boundary. Then, when the algorithm of interest sweeps over the interior matrix, it can access nodes that are one or two steps away in _any_ direction from all the interior nodes. This simplifies coding as the alternative would be separate loops for the sides with modified algorithms that builds the boundary conditions into the operations. In GMT, we simply have one function that can be called to set the padding values from the boundary condition settings (should there not be data available). Note that the padding, while typically 2, can be changed internally, and there are places where a pad of 1 is used.

**What are the boundary conditions?**

GMT supports both natural, periodic, geographic, and data boundary conditions. _Natural_ means the curvature at the edge is zero, while _periodic_ means we have a repeating boundary (360 longitude wrapping is such a boundary condition, where _lon = 361_ is same as _lon = 1_). _Geographic_ is a special mix of periodic in longitude but pole-wrapping in latitude (as we go over the pole along a meridian, longitude changes by 180). Finally, _data_ boundary condition means we are extracting a subset of a data grid, but since it is a subset (and hence data beyond it are available) we can use real data to fill in the pad instead of applying some mathematical condition. Now, when a derivative is required at the boundary it will be data-based. This is clearly superior to _any other scheme_.

**OK, make sense, so what is the problem?**

Once the basic input/output of grids in GMT are not files but memory references to grids maintained by Matlab, Python, or Julia, then two situations arise that require attention. **Input:** Receiving a grid from an external environment is the equivalent of reading a grid file, except now there may be _no padding_ associated with the grid. If the module that receives this pad-less grid actually requires padding for its algorithm to work, and the input is read-only memory, we must _duplicate_ the grid and place it in a new padded matrix and do the work on that padded matrix. This increases memory requirements. **Output:** At the point of writing the output, the internal matrix is typically padded, but the receiver (say, Julia or Python) has no need for a pad nor expects it. Thus, sending padded output from GMT to externals is generally not appreciated and leads to coding complications for the external developer.

**Which GMT modules require a grid pad?**

While some modules that create a grid from other input (e.g., **surface**) needs padding, here I am only concerned with the modules that take grid input, since that is where the padding issue will be present. These include **grdimage** and **grdview** (for projections to properly fill new nodes near boundaries), **grdproject** (same reason), **grdsample** (2-D bicubic sampling ), **grdtrack** (involves 2-D bicubic sampling). I am not sure if **grdcontour** may also be in this category, and there may be others [forcing _GMT_Create_Session_ to set pad to 0, recompile, and run tests might tell us more!]. Also, note that some GMT modules call other modules, for instance **grdblend** may call **grdsample**, and thus inherits the need for a pad. This does not change the conclusions that (a) many grid-reading modules do not require boundary conditions (e.g., **grdinfo**, **grdcut**, **grdclip**, etc.) and (b) those that do really require it and rewriting the code to _not_ work on matrices with pads would be extensive, bug-prone, and not likely a good investment of time.

**So how are things working now with memory grids coming in?**

While we are still making adjustments to deal with this relatively new scenario, by and large GMT itself is coping with the input. Those modules that must have padding will now check if input has a pad and if not we do the data duplication/expansion, work on the temporary matrix instead, and free it when done with it. This is largely unnoticeable by the user as we are not talking huge amounts of memory unless the input is truly large and the user's environment is memory-limited.

**How are things with externals receiving GMT-produced grids by memory?**

I suspect this is where there is room for improvement. Regardless of the need to use padding internally to compute results, there seems to be no good reason to return a _padded_ array back to the calling environment. Thus, when _GMT_Write_Data_ is called and the output is a grid destined fo a memory location, the pad should always be stripped off. If that was true then the external user would never notice the padding issue.

**Are there downsides of not having pads outside GMT?**

Yes there are. Consider an external tool that extracts the subregion 0/20/0/20 from a global topography grid, then passes that in memory via a call to **grdtrack** to sample it along lines that crosses the domain. Because of no padding, **grdtrack** has to apply natural boundary conditions during the interpolation along tracks. The interpolated values very close to the grid boundary will generally differ from what the GMT command-line interface (CLI) equivalent reading that subset would do since the pad would have actual data, not pseudo-values derived from the natural boundary condition. So the CLI result and the external call to **grdtrack** would differ in subtle ways, with the CLI being the correct result. Of course, an informed user could compensate and ask for -1/22/-2/22 as the subset, do the interpolation, and then crop to the desired window, but this involves more work as well as considerable awareness of the effects we have discussed. With the GMT developers historically being very concerned with "doing the right thing" and "screw people with low-RAM computers", you can see why we rather want an exact result than an almost exact one.

Hopefully the @GenericMappingTools/core can comment on this, and if something is incorrect or missing please let me know so I can update this description.