Maybe not an answer you were expecting, it appears PDL does not provide interface to some GSL functions. The PDL::GSL::INTERP describes 1D interpolation only. However, there are 2D interpolation functions for irregular/non-uniform X-Y grid, exactly as you require. Then, after interpolating for points on regular grid, you could proceed to visualization on a 2D image, etc.
There is "gsl_interp2d.h" in my Strawberry Perl installation, so it must be possible to call these functions from existing libgsl dll. Writing PDL bindings is beyond my scope, but perhaps it's not prohibitively difficult to use FFI in this case.
Another possibility, I had in mind even before finding the above, is that irregular grid is just a specific case of any random set (not a grid) of (X,Y) points for which f(X,Y) is known. From f(X,Y) for any set, to f on regular grid, for reasonably well-behaving function and if linear interpolation is good enough -- I did that before, though it's not PDL. First, triangulate X-Y plane, then interpolate for regular grid points inside each triangle.