Use map instead of array in C++ to protect searching outside of array bounds?

My advice is to replace all read accesses to the array by a function. For example, arr[i][j] by getarr(i,j). That way, all your algorithmic code stays more or less unchanged and you can easily return NODATA for indices outside bounds.

But I must admit that it is only my opinion.

I've had to do this before and the fastest solution was to expand the region with NODATA values and iterate over the interior. This way the core loop is simple for the compiler to optimize.

If this is not a computational hot-spot in the code, I'd go with Serge's approach instead though.

To minimize rippling effects I used an array structure with explicit row/column strides, something like this:

class Grid {
  private:
    shared_ptr<vector<double>> data;
    int origin;
    int xStride;
    int yStride;

  public:
    Grid(int nx, int ny) : 
      data( new vector<double>(nx*ny) ),
      origin(0),
      xStride(1), 
      yStride(nx) {
    }

    Grid(int nx, int ny, int padx, int pady) :
      data( new vector<double>((nx+2*padx)*(ny+2*pady));
      xStride(1), 
      yStride(nx+2*padx),
      origin(nx+3*padx) {
    }

  double& operator()(int x, int y) {
     return (*data)[origin + x*xStride + y*yStride];
  }
}

Now you can do

  Grid g(5,5,1,1);
  Grid g2(5,5);
  //Initialise
  for(int i=0; i<5; ++i) {
    for(int j=0; j<5; ++j) {
     g(i,j)=i+j;
    }
  }
  // Convolve (note we don't care about going outside the
  // range, and our indices are unchanged between the two
  // grids.
  for(int i=0; i<5; ++i) {
    for(int j=0; j<5; ++j) {
      g2(i,j)=0;
      g2(i,j)+=g(i-1,j);
      g2(i,j)+=g(i+1,j);
      g2(i,j)+=g(i,j-1);
      g2(i,j)+=g(i,j+1);
    }
  }

Aside: This data structure is awesome for working with transposes, and sub-matrices. Each of those is just an adjustment of the offset and stride values.

Solution 1 is the standard solution. It takes maximum advantage of modern computer architectures, where a few bytes of memory are no big deal, and correct instruction prediction accelerates performance. As you keep accessing memory in a predictable pattern (with fixed strides), the CPU prefetcher will successfully read ahead.

Solution 2 saves a small amount of memory, but the special handling of the edges incurs a real slowdown. Still, the large chunk in the middle benefits from the prefetcher.

Solution 3 is horrible. Map access is O(log N) instead of O(1), and in practice it can be 10-20 times slower. Maps have poor locality of reference; the CPU prefetcher will not kick in.

If simple means "easy to read" I'd recommend you declare a class with an overloaded [] operator. Use it like a regular array but it'll have bounds checking to handle NODATA.

If simple means "high performance" and you have sparse grid with isolated DATA consider implementing linked lists to the DATA values and implement optimal operators that go directly to tge DATA values.

1 wastes memory proportional to your overall rectangle size, 3/maps are clumsy here, 2 is actually very easy to do:

T d[X][Y] = ...;

for (int x = 0; x < X; ++x)
    for (int y = 0; y < Y; ++y)  // move over d[x][y] centres
    {
        T r[3][3] = { { d[i,j], d[i,j], d[i,j] },
                        d[i,j], d[i,j], d[i,j] },
                        d[i,j], d[i,j], d[i,j] } };

        for (int i = std::min(0, x-1); i < std::max(X-1, x+1); ++i)
            for (int j = std::min(0, y-1); j < std::max(Y-1, y+1); ++j)
                if (d[i][j] != NoData)
                    r[i-x][j-y] = d[i][j];

        // use r for whatever...
    }

Note that I'm using signed int very deliberately so x-1 and y-1 don't become huge positive numbers (as they would with say size_t) and break the std::min logic... but you could express it differently if you had some reason to prefer size_t (e.g. x == 0 ? 0 : x - 1).