This is not supposed to be a solution, but rather an explanation and a possible, although ugly workaround.
From time to time there are reports to the gnuplot mailinglists about this issue, but it seems to be related to the viewers. It has to do with the way gnuplot creates the surfaces plots. These are drawn as polygons, which are stitched together. The Moiré patterns you are showing come from wrong rendering between two polygons. That depends on the viewer, the viewer settings and the zoom factor.
The easiest example, to show that effect is the following Postscript file:
%!PS-Adobe-2.0
50 50 moveto 50 0 rlineto 0 50 rlineto -50 0 rlineto closepath 0 setgray fill
100 50 moveto 50 0 rlineto 0 50 rlineto -50 0 rlineto closepath 0 setgray fill
Save this file e.g. as moire.ps and view it, or convert it with ps2pdf and view it. With the Acrobat reader 9.5.1 I see the following:
The Acrobat Reader has a setting Preferences -> Page Display -> Enhance thin lines which can prevent this problem, but causes problems on other parts.
On my system (Debian), all viewers show this patterns, mupdf, firefox, ghostscript, pdftocairo, libpoppler` etc.
So, what to do? For myself I use the following workaround. I splot to a png with high resolution, and reread that file later with plot ... with rgbimage. Then you get your heatmap as bitmap, and the rest is vectorial. In most cases this is no problem, because in any way you have some measurement data with limited resolution, which you interpolate.
Based on the question gnuplot contour line color: set style line and set linetype not working, here is how you can implement it:
reset
set lmargin at screen 0.05
set rmargin at screen 0.85
set bmargin at screen 0.1
set tmargin at screen 0.9
set pm3d map interpolate 20,20
unset key
set cntrparam bspline
set cntrparam points 10
set cntrparam levels increment -6,-6,-24
set contour surface
set linetype 1 lc rgb "blue" lw 2
set linetype 2 lc rgb "blue"
set linetype 3 lc rgb "black"
set linetype 4 lc rgb "orange"
set linetype 5 lc rgb "yellow"
set palette rgb 33,13,10 #rainbow (blue-green-yellow-red)
set cbrange [-18:0]
unset border
unset xtics
unset ytics
set angles degree
r = 3.31
set xrange[-r:r]
set yrange[-r:r]
set colorbox user origin 0.9,0.1 size 0.03,0.8
##################### start changes ##############
set autoscale fix
RES_X = 2000
RES_Y = 2000
save('settings.tmp')
set lmargin at screen 0
set rmargin at screen 1
set bmargin at screen 0
set tmargin at screen 1
unset colorbox
set terminal pngcairo size RES_X, RES_Y
set output '3d-polar-inc.png'
splot 'new_test.dat' nocontour
unset output
load('settings.tmp')
# mapping of the coordinates for the png plotting later
X0 = GPVAL_X_MIN
Y0 = GPVAL_Y_MIN
DX = (GPVAL_X_MAX - GPVAL_X_MIN)/real(RES_X)
DY = (GPVAL_Y_MAX - GPVAL_Y_MIN)/real(RES_Y)
C0 = GPVAL_CB_MIN
DC = GPVAL_CB_MAX - GPVAL_CB_MIN
C(x) = (x/255.0) * DC + C0
# now plot the png
#set terminal pdfcairo size 10cm,10cm
#set output '3d-polar.pdf'
set terminal postscript eps color level3 size 10cm,10cm solid
set output '3d-polar-eps.eps'
set multiplot
set cbrange[GPVAL_CB_MIN:GPVAL_CB_MAX]
plot '3d-polar-inc.png' binary filetype=png \
origin=(X0, Y0) dx=DX dy=DY \
using (C($1)):(C($2)):(C($3)) \
with rgbimage, \
NaN with image t '' # hack for getting the colorbox
# plot the contours
unset surface
unset pm3d
splot 'new_test.dat' w l
###################### end changes #################
# now plot the polar grid only
set style line 11 lc rgb 'black' lw 2 lt 0
set grid polar ls 11
set polar
set logscale r 10
set rrange[10:20000]
unset raxis
set rtics format '' scale 0
#set rtics axis scale
set rtics (20,50,100,200,500,1000,2000,5000,10000,20000)
do for [i=-150:180:30] {
dum = r+0.15+0.05*int(abs(i/100))+0.05*int(abs(i/140))-0.05/abs(i+1)
set label i/30+6 at first dum*cos(i), first dum*sin(i) center sprintf('%d', i)
}
set label 20 at first 0, first -(log(20)/log(10)-1) center "20"
set label 100 at first 0, first -(log(100)/log(10)-1) center "100"
set label 200 at first 0, first -(log(200)/log(10)-1) center "200"
set label 1000 at first 0, first -(log(1000)/log(10)-1) center "1k"
set label 2000 at first 0, first -(log(2000)/log(10)-1) center "2k"
set label 10000 at first 0, first -(log(10000)/log(10)-1) center "10k"
set label 20000 at first 0, first -(log(20000)/log(10)-1) center "20k"
plot NaN w l
unset multiplot
unset output
With pdfcairo this gives a 1.7 MB pdf file, with epslatex level3 (this option is available in the 4.7 development version only) you get a 1.5 MB eps file, which can be converted with epstopdf to a 136 KB pdf file.
See also my answer to Big data surface plots: Call gnuplot from tikz to generate bitmap and include automatically? on TeX.SX.