gnuplot contour plot hatched lines

Question

I'm using gnuplot for contour plot of a several function. This is for optimization problem. I have 3 functions:

1. f(x,y)
2. g1(x,y)
3. g2(x,y)

both g1(x,y) and g2(x,y) are constraints and would like to plot on top of the contour plot of f(x,y).

Here is the textbook example:

Here is my attempt to replicate it in gnuplot, thanks to @theozh.

### contour lines with labels
reset session

f(x,y)=(x**2+y-11)**2+(x+y**2-7)**2
g1(x,y)=(x-5)**2+y**2
g2(x,y) = 4*x+y

set xrange [0:6]
set yrange [0:6]
set isosample 250, 250
set key outside

set contour base
set cntrparam levels disc 10,30,75,150,300,500,850,1500 
unset surface
set table $Contourf
    splot f(x,y)
unset table

set contour base
set cntrparam levels disc 26
unset surface
set table $Contourg1
    splot g1(x,y)
unset table

set contour base
set cntrparam levels disc 20
unset surface
set table $Contourg2
    splot g2(x,y)
unset table

set style textbox opaque noborder

set datafile commentschar " "
plot for [i=1:8] $Contourf u 1:2:(i) skip 5 index i-1 w l lw 1.5 lc var title columnheader(5)
replot $Contourg1 u 1:2:(1) skip 5 index 0 w l lw 4 lc 0 title columnheader(5)
replot $Contourg2 u 1:2:(1) skip 5 index 0 w l lw 4 lc 0 title columnheader(5)

I would like to replicate the textbook picture in the gnuplot example. How to do a hatch mark on the functions g1 and g2, the thick black line in plot above.

@theozh provided an excellent solution below. However, the method doesnot work for steep curves. As an example

reset session
unset key

set size square

g(x,y) = -0.8-1/x**3+y

set xrange [0:4]
set yrange [0:4]
set isosample 250, 250
set key off

set contour base
unset surface

set cntrparam levels disc 0
set table $Contourg
    splot g(x,y)
unset table

set angle degree
set datafile commentschar " "

plot $Contourg u 1:2 skip 5 index 0 w l lw 2 lc 0 title columnheader(5)

set style fill transparent pattern 4
replot $Contourg u 1:2:($2+0.2) skip 5 index 0 w filledcurves lc 0 notitle 

yields the following figure. Is there a way to use different offsets, for example offset x values for x < 1.3 and for x > 1.3 offset y values. This would yield a much better filled curve. A matlab implementations of what I was looking for can be found here: https://www.mathworks.com/matlabcentral/fileexchange/29121-hatched-lines-and-contours.

In replcating @Ethans program, I get the following, the dashtype is relatively thick compared to @Ethan not sure why, I'm using gnuplot v5.2 and wxt terminal.

When I replicate @theozh code, it works very well except for closed contours, not sure why? see below for example:

f(x,y)=x*exp(-x**2-y**2)+(x**2+y**2)/20
g1(x,y)= x*y/2+(x+2)**2+(y-2)**2/2-2

set xrange [-7:7]
set yrange [-7:7]
set isosample 250, 250
set key outside

set contour base
unset surface

set cntrparam levels disc 4,3.5,3,2.5,2,1.5,1,0.5,0 
set table $Contourf
    splot f(x,y)
unset table

set cntrparam levels disc 0
set table $Contourg1
    splot g1(x,y)
unset table

# create some extra offset contour lines
# macro for setting contour lines
ContourCreate = '\
    set cntrparam levels disc Level; \
    set table @Output; \
        splot @Input; \
    unset table'

Level = 0.45
Input = 'g1(x,y)'
Output = '$Contourg1_ext'
@ContourCreate


# Macro for ordering the datapoints of the contour lines which might be split
ContourOrder = '\
stats @DataIn skip 6 nooutput; \
N = STATS_blank-1; \
set table @DataOut; \
    do for [i=N:0:-1] { plot @DataIn u 1:2 skip 5 index 0 every :::i::i with table }; \
unset table'

DataIn = '$Contourg1'
DataOut = '$Contourg1_ord'
@ContourOrder

DataIn = '$Contourg1_ext'
DataOut = '$Contourg1_extord'
@ContourOrder


# Macro for reversing a datablock
ContourReverse = '\
set print @DataOut; \
    do for [i=|@DataIn|:1:-1] { print @DataIn[i]}; \
set print'

DataIn = '$Contourg1_extord'
DataOut = '$Contourg1_extordrev'
@ContourReverse

# Macro for adding datablocks
ContourAdd = '\
set print @DataOut; \
    do for [i=|@DataIn1|:1:-1] { print @DataIn1[i]}; \
    do for [i=|@DataIn2|:1:-1] { print @DataIn2[i]}; \
set print'

DataIn1 = '$Contourg1_ord'
DataIn2 = '$Contourg1_extordrev'
DataOut = '$Contourg1_add'
@ContourAdd


set style fill noborder 
set datafile commentschar " "
plot \
    for [i=1:8] $Contourf u 1:2:(i) skip 5 index i-1 w l lw 1.5 lc var title columnheader(5), \
    $Contourg1 u 1:2 skip 5 index 0 w l lw 2 lc 0 title columnheader(5), \
    $Contourg1_add u 1:2 w filledcurves fs transparent pattern 5 lc rgb "black" notitle

Answer 1

Another possibility is to use a custom dash pattern, as shown below: By the way, it is almost never correct to use "replot" to compose a single figure.

# Additional contour levels displaced by 0.2 from the original
set contour base
set cntrparam levels disc 20.2
unset surface
set table $Contourg2d
    splot g2(x,y)
unset table

set contour base
set contour base
set cntrparam levels disc 26.2
unset surface
set table $Contourg1d
    splot g1(x,y)
unset table

set linetype 101 lc "black" linewidth 5 dashtype (0.5,5)

plot for [i=1:8] $Contourf u 1:2:(i) skip 5 index i-1 w l lw 1.5 lc var title columnheader(5), \
        $Contourg1 u 1:2:(1) skip 5 index 0 w l lw 1 lc "black" title columnheader(5), \
        $Contourg2 u 1:2:(1) skip 5 index 0 w l lw 1 lc "black" title columnheader(5), \
        $Contourg1d u 1:2:(1) skip 5 index 0 w l linetype 101 notitle, \
        $Contourg2d u 1:2:(1) skip 5 index 0 w l linetype 101 notitle

Amended to show use of contours offset so that the dashes are only on one side of the line.

Answer 2

Here is the solution you (and I) were hoping for. You just enter the hatch parameters into a datablock: TiltAngle in degrees (>0°: left side, <0° right side in direction of path), HatchLength and HatchGap in pixels. The procedure has become a bit lengthy but it does what you want. I have tested it with gnuplot 5.2.8 and 5.4.1 and wxt and qt terminal.

What the procedure basically does:

1. determines the angle between two consecutive points of the data input curve
2. interpolates datapoints along the curve according to HatchSeparation
3. Scales everything such that is independent of graph scale and terminal size (this, however, requires a dummy plot of the curves without hatch lines for getting the gnuplot variables GPVAL_X_MAX, GPVAL_X_MIN, GPVAL_TERM_XMAX, GPVAL_TERM_XMIN, GPVAL_Y_MAX, GPVAL_Y_MIN, GPVAL_TERM_YMAX, GPVAL_TERM_YMIN.

Limitations:

• does not work (yet) with logarithmic axes
• several independent paths need to be separated by two empty lines
• if you are using it together with your contour lines you have to make sure that the contour line datapoints are in the right order (see comment in my first answer). Furthermore, gnuplot might generate contour lines which are separated with a single empty line. I can address this if I find some time or/and somebody requests it.

Edit: (completely revised version)

The previous script (to my opinion) was pretty messy and difficult to follow (although nobody complained ;-). I removed the calls to subprocedures and hence the prefixes for variables in the subprocedures and put all in one script, except the test data generation.

Have fun with hatching your lines! Comments and improvements are welcome!

Test data generation: SO57118566_createTestData.gp

### Create some circle test data

FILE = "SO57118566.dat"
set angle degrees

# create some test data
# x     y     r     a0  a1    N
$myCircleParams <<EOD
 1.0    0.3   0.6   0   360   120
 2.4    0.3   0.6   0   360   120
 3.8    0.3   0.6   0   360   120
 1.7   -0.3   0.6   0   360   120
 3.1   -0.3   0.6   0   360   120
EOD

X(n)  = real(word($myCircleParams[n],1))   # center x
Y(n)  = real(word($myCircleParams[n],2))   # center y
R(n)  = real(word($myCircleParams[n],3))   # radius
A0(n) = real(word($myCircleParams[n],4))   # start angle
A1(n) = real(word($myCircleParams[n],5))   # end   angle
N(n)  =  int(word($myCircleParams[n],6))   # number of samples

set table FILE
    do for [i=1:|$myCircleParams|] {
        set samples N(i)
        plot [A0(i):A1(i)] '+' u (X(i)+R(i)*cos($1)):(Y(i)+R(i)*sin($1))
    }
unset table

set size ratio -1
plot FILE u 1:2:-2 w l lc var
### end of script

Strange enough, the previous version worked for gnuplot5.2.0 to 5.2.7 but not for gnuplot>=5.2.8. With this current script it is vice versa, but I haven't yet found out why.

Update: Finally found why it wasn't working with <=5.2.7. Apparently something with the scaling which has changed between 5.2.7 and 5.2.8. Other terminals than wxt or qt might have different scaling factors. You need to add/change the lines (already added in the script below):

Factor = GPVAL_VERSION==5.2 && int(GPVAL_PATCHLEVEL)<=7 ? \
         GPVAL_TERM eq "wxt" ? 20 : GPVAL_TERM eq "qt" ? 10 : 1 : 1
Rxaupu = (GPVAL_X_MAX-xmin)/(GPVAL_TERM_XMAX-xtmin)*Factor   # x ratio axes units to pixel units
Ryaupu = (GPVAL_Y_MAX-ymin)/(GPVAL_TERM_YMAX-ytmin)*Factor   # y

Script: (tested with gnuplot 5.2.0, 5.2.7, 5.2.8, 5.4.1)

### Add hatch pattern to a curve
reset session

FILE = "SO57118566.dat"

set size ratio -1       # set same x,y scaling
set angle degree
unset key

# plot path without hatch lines to get the proper gnuplot variables: GPVAL_...
plot FILE u 1:2:-2 w l lc var

# Hatch parameters:
# TiltAngle     >0°: left side, <0° right side in direction of path
# HatchLength   hatch line length in pixels
# HatchGap      separation of hatch lines in pixels
# TA   HL   HG   Color
$myHatchParams <<EOD
 -90   10    5   0x0080ff
 -30   15   10   0x000000
  90    5    3   0xff0000
  45   25   12   0xffff00
 -60   10    7   0x00c000
EOD
# extract hatch parameters
TA(n)    = real(word($myHatchParams[n],1))   # TiltAngle
HL(n)    = real(word($myHatchParams[n],2))   # HatchLength
Gpx(n)   = real(word($myHatchParams[n],3))   # HatchGap in pixels
Color(n) =  int(word($myHatchParams[n],4))   # Color

# terminal constants
xmin   = GPVAL_X_MIN
ymin   = GPVAL_Y_MIN
xtmin  = GPVAL_TERM_XMIN
ytmin  = GPVAL_TERM_YMIN
Factor = GPVAL_VERSION==5.2 && int(GPVAL_PATCHLEVEL)<=7 ? \
         GPVAL_TERM eq "wxt" ? 20 : GPVAL_TERM eq "qt" ? 10 : 1 : 1
Rxaupu = (GPVAL_X_MAX-xmin)/(GPVAL_TERM_XMAX-xtmin)*Factor   # x ratio axes units to pixel units
Ryaupu = (GPVAL_Y_MAX-ymin)/(GPVAL_TERM_YMAX-ytmin)*Factor   # y

Angle(dx,dy) = dx==0 && dy==0 ? NaN : atan2(dy,dx)    # -180° to +180°, NaN if dx,dy==0
LP(dx,dy)    = sqrt(dx**2 + dy**2)        # length of path segment
ax2px(x)     = (x-xmin)/Rxaupu  + xtmin   # x axes coordinates to pixel coordinates
ay2py(y)     = (y-ymin)/Rxaupu  + ytmin   # y
px2ax(x)     = (x-xtmin)*Rxaupu + xmin    # x pixel coordinates to axes coordinates
py2ay(y)     = (y-ytmin)*Rxaupu + ymin    # y

# create datablock $Path with pixel coordinates and cumulated path length
stats FILE u 0 nooutput    # get number of blocks of input file
N = STATS_blocks
set table $Path
    do for [i=0:N-1] {
        x1 = y1 = NaN
        Length = 0
        plot FILE u (x0=x1, x1=ax2px($1)):(y0=y1, y1=ay2py($2)): \
            (dx=x1-x0, dy=y1-y0, ($0>0?Length=Length+LP(dx,dy):Length)) index i w table
        plot '+' u ('') every ::0::1 w table    # two empty lines
    }
unset table

# create hatch lines table
# resample data in equidistant steps along the length of the path
$Temp <<EOD    # datablock $Temp definition required for function definition below
EOD
x0(n)    = real(word($Temp[n],1))                # x coordinate
y0(n)    = real(word($Temp[n],2))                # y coordinate
r0(n)    = real(word($Temp[n],3))                # cumulated path length 
ap(n)    = Angle(x0(n+1)-x0(n),y0(n+1)-y0(n))    # path angle
ah(n,i)  = ap(n)+TA(i+1)                         # hatch line angle
Frac(n)  = (ri-r0(n))/(r0(n+1)-r0(n))            # interpolation along
hsx(n)   = (x0(n) + Frac(n)*(x0(n+1)-x0(n)))     # x hatch line start point
hsy(n)   = (y0(n) + Frac(n)*(y0(n+1)-y0(n)))     # y
hex(n,i) = (hsx(n) + HL(i+1)*cos(ah(n,i)))       # x hatch line end point
hey(n,i) = (hsy(n) + HL(i+1)*sin(ah(n,i)))       # y

# create datblock with hatchlines x,y,dx,dy
set print $HatchLines
    do for [i=0:N-1] {
        set table $Temp
            splot $Path u 1:2:3 index i
        unset table

        ri = -Gpx(i+1)
        do for [j=1:|$Temp|-2] {
            if (strlen($Temp[j])==0 || $Temp[j][1:1] eq '#') {print $Temp[j]}
            else {
                while (ri<r0(j)) {
                    ri = ri + Gpx(i+1)
                    print sprintf("%g %g %g %g", \
                    xs=px2ax(hsx(j)), ys=py2ay(hsy(j)), \
                    px2ax(hex(j,i))-xs, py2ay(hey(j,i))-ys)
                }
            }
        }
        print ""; print ""   # two empty lines
    }
set print

plot $Path u (px2ax($1)):(py2ay($2)):(Color(column(-2)+1))  w l lc rgb var, \
     $HatchLines u 1:2:3:4:(Color(column(-2)+1)) w vec nohead lc rgb var
### end of script

Result: