Although a trivial(?) task for a seasoned gmt user, I think this would be a very nice addition to grdcut. It is a fairly common operation. Saving time is always nice.
One thing though; retaining all values inside the polygon and set the outside to NaN is nice, but shouldnt there also be calculated a new -R based on the polygon? I can see many cases where using the original -R from the input grid, would cause a too big region to be used, containing only NaN’s.