Trying to run the @theozh script How to display non-rectangular and ungridded data as a map? I've gotten the very same plot except for straight lines that appear enclosing the various colored parts. I have run the same script (see below) except for the first lines since Windows CoreUtils sort command doesn't work well in my case, so I have generated the sorted $Data file otherwise. But from "# definition of quite a few variables..." on, the script is identical. Does anyone know how to get rid of those enclosing lines? Maybe$Datashould include blank lines between rows ? (that might be the only difference).
### Delaunay triangulation (gnuplot implementation attempt, requires gnuplot 5.4) by theozh
reset session
# get some test data
Random=0 # set to 0 to read data from existing FILE
if (Random) {
FILE = "tbTriangulationRandom.dat"
set print FILE
do for [i=0:99] {
print sprintf("%g %g %g",x=invnorm(rand(0))*10,y=invnorm(rand(0))*10,x*y)
}
set print
}
else {
FILE = "tbTriangulationTest.dat" # IT COMES IN THE LINK PROVIDED BEFORE
}
# sort data by x increasing values and if x is identical by increasing y values
set format y ""
set table $Data0
plot FILE u 1:2:1 smooth zsort
set table $Data3
plot FILE u 2:2:1 smooth zsort
set table $Data33
plot FILE u 3:2:1 smooth zsort
unset table
set print $Data4
do for [i=1:|$Data0|] {
print $Data0[i],$Data3[i],$Data33[i]
if (i<|$Data0|) { if (word($Data0[i],1) ne word($Data0[i+1],1)) { print "" } }
}
set print
set table $Data0
plot $Data4 u 1:3:3 smooth zsort
set table $Data3
plot $Data4 u 3:3:3 smooth zsort
set table $Data33
plot $Data4 u 5:3:3 smooth zsort
unset table
set print $Data44
do for [i=1:|$Data0|] { if ($Data4[i] ne "") { print $Data0[i],$Data3[i],$Data33[i]; print " "
}}
set print
set table $Data
plot $Data44 u 1:3:5 w table
unset table
set format y
# definition of quite a few variables and functions
colX = 1 # x column
colY = 2 # y column
colZ = 3 # z column
N = |$Data| # number of points
EDGES = '' # list of (inner) edges
HULL = '' # list of hull segments
TRIANGLES = '' # list of triangles
HULLPOINTS = '' # list of hullpoints
array Px[N] # set point array size
array Py[N] # set point array size
array Pz[N] # set point array size
newEdge(p1,p2) = sprintf(" %d %d ",p1,p2)
Edge(n) = sprintf(" %s %s ",word(EDGES,2*n-1),word(EDGES,2*n))
EdgeP(n,p) = int(word(EDGES,2*n-2+p))
changeEdge(n,p1,p2) = (_edge=Edge(n), _pos = strstrt(EDGES,_edge), _pos ? \
EDGES[1:_pos-1].newEdge(p1,p2). \
EDGES[_pos+strlen(_edge):strlen(EDGES)] : EDGES)
TriangleCount(n) = words(TRIANGLES)/3
TriangleN(n) = sprintf(" %s %s %s ", \
word(TRIANGLES,3*n-2),word(TRIANGLES,3*n-1),word(TRIANGLES,3*n))
newTriangle(p1,p2,p3) = p1<p2 && p2<p3 ? sprintf(" %d %d %d ",p1,p2,p3) : \
p1<p3 && p3<p2 ? sprintf(" %d %d %d ",p1,p3,p2) : \
p2<p1 && p1<p3 ? sprintf(" %d %d %d ",p2,p1,p3) : \
p2<p3 && p3<p1 ? sprintf(" %d %d %d ",p2,p3,p1) : \
p3<p1 && p1<p2 ? sprintf(" %d %d %d ",p3,p1,p2) : \
sprintf(" %d %d %d ",p3,p2,p1)
changeTA(n,p1,p2,p3) = (TA=TriangleN(n), _pos = strstrt(TRIANGLES,TA), _pos ? \
TRIANGLES[1:_pos-1].newTriangle(p1,p2,p3). \
TRIANGLES[_pos+strlen(TA):strlen(TRIANGLES)] : TRIANGLES)
TAp(n,p) = int(word(TRIANGLES,3*n-3+p))
TAx(n,p) = Px[TAp(n,p)] # x-coordinate of point p of triangle n
TAy(n,p) = Py[TAp(n,p)] # y-coordinate of point p of triangle n
HullP(n,p) = int(word(HULL,2*n-2+p)) # hull segment point number
HScount(n) = int(words(HULL))/2 # number of hull segments
getHullPoints(n) = (_tmp = '', sum [_i=1:words(HULL)] ((_s=' '.word(HULL,_i).' ', _tmp = \
strstrt(_tmp,_s) ? _tmp : _tmp._s ), 0), _tmp)
removeFromHull(seg) = (seg, _pos = strstrt(HULL,seg), _pos ? \
HULL[1:_pos-1].HULL[_pos+strlen(seg):strlen(HULL)] : HULL)
# orientation of 3 points, either -1=clockwise, 0=collinear, 1=counterclockwise
Orientation(p1,p2,p3) = sgn((Px[p2]-Px[p1])*(Py[p3]-Py[p1]) - (Px[p3]-Px[p1])*(Py[p2]-Py[p1]))
# check for intersection of segment p1-p2 with segment p3-p4, 0=no intersection, 1=intersection
IntersectCheck(p1,p2,p3,p4) = (Orientation(p1,p3,p2)==Orientation(p1,p4,p2)) || \
(Orientation(p3,p1,p4)==Orientation(p3,p2,p4)) ? 0 : 1
Sinus(p1,p2) = (_dx=Px[p2]-Px[p1], _dy=Py[p2]-Py[p1], _dy/sqrt(_dx**2 + _dy**2))
### Macros for later use
# Creating inner edges datablock macro
CreateEdgeDatablock = '\
set print $EDGES; \
do for [i=1:words(EDGES)/2] { \
p1 = int(word(EDGES,2*i-1)); \
p2 = int(word(EDGES,2*i)); \
print sprintf("%g %g %g %g %d %d",Px[p1],Py[p1],Px[p2]-Px[p1],Py[p2]-Py[p1],p1,p2) \
}; \
set print '
# Creating hull datablock macro
CreateHullDatablock = '\
set print $HULL; \
do for [i=1:words(HULL)/2] { \
p1 = int(word(HULL,2*i-1)); \
p2 = int(word(HULL,2*i)); \
print sprintf("%g %g %g %g %d %d",Px[p1],Py[p1],Px[p2]-Px[p1],Py[p2]-Py[p1],p1,p2) \
}; \
set print '
# plotting everything
PlotEverything = '\
plot $EDGES u 1:2:3:4 w vec lw 1.0 lc "red" nohead, \
$HULL u 1:2:3:4 w vec lw 1.5 lc "blue" nohead, \
$Data u 1:2 w p pt 7 ps 0.5 lc "black", \
$Data u 1:2:($0+1) w labels offset 0.5,0.5 '
# put datpoints into arrays
set table $Dummy
plot $Data u (Px[int($0)+1]=column(colX),Py[int($0)+1]=column(colY),Pz[int($0)+1]=column(colZ),'') w table
unset table
# get first m>=3 points which are not all collinear
HULL = HULL.newEdge(1,2) # add first 2 points to hull in any case
do for [p=3:N] {
if (Orientation(p-2,p-1,p)==0) { # orientation==0 if collinear
HULL = HULL.newEdge(p-1,p)
}
else { break } # stop if first >=3 non-collinear points found
}
HPcountInit = words(getHullPoints(0)) # get initial number of hull points
# actual plotting starts from here
set offset 1,1,1,1
set key noautotitle
set palette rgb 33,13,10
set rmargin screen 0.8
plot $Data u 1:2 w p pt 7 ps 0.5 lc "black", \
'' u 1:2:($0+1) w labels offset 0.5,0.5
set label 1 at graph 0.02,0.97 "Adding points... "
# loop all data points
do for [p=HPcountInit+1:N] {
print sprintf("### Adding P%d",p)
HPlist = getHullPoints(0)
HPcount = words(HPlist)
set print $NewConnections # initalize/empty datablock for new connections
print ""
set print
do for [hpt in HPlist] { # loop and check all hull points
hp = int(hpt)
# print sprintf("Check hull point P%d", hp)
c = 0
do for [hs=1:HScount(0)] { # loop all hull segments
hp1 = HullP(hs,1)
hp2 = HullP(hs,2)
# print sprintf("Check %d-%d with %d-%d", hp1, hp2, hp, p)
if (p!=hp1 && p!=hp2 && hp!=hp1 && hp!=hp2) {
c = c || IntersectCheck(hp1,hp2,hp,p)
if (c) { break }
}
}
if (c==0) { # if no intersections with hull
set print $NewConnections append # add new connections to datablock
print sprintf("%g %g", hp, Sinus(p,hp))
set print
}
}
# sort datablock clockwise (a bit cumbersome in gnuplot)
set table $ConnectSorted
plot $NewConnections u 1:2:2 smooth zsort # requires gnuplot 5.4.0
set table $Dummy
plot Connect='' $ConnectSorted u (Connect=Connect.sprintf(" %d",$1),'') w table
unset table
# add new edges
Ccount = int(words(Connect))
do for [i=1:Ccount-1] {
p1 = int(word(Connect,i))
p2 = int(word(Connect,i+1))
TRIANGLES = TRIANGLES.sprintf(" %d %d %d ", p1<p2?p1:p2, p2<p1?p1:p2, p) # numbers in ascending order
if (i==1) { HULL=HULL.newEdge(p1,p) }
if (i>1 && i<Ccount) { EDGES = EDGES.newEdge(p1,p) }
if (i==Ccount-1) {
HULL = HULL.newEdge(p2,p)
}
if (p!=HPcountInit+1) { # remove hull segments, except initial ones
NewEdge = p1<p2 ? sprintf(" %d %d ",p1,p2) : sprintf(" %d %d ",p2,p1)
# print sprintf("Remove %s",NewEdge)
HULL = removeFromHull(NewEdge)
EDGES = EDGES.NewEdge
@CreateEdgeDatablock
@CreateHullDatablock
@PlotEverything
}
}
}
# flip diagonal of a quadrangle if Det(p1,p2,p3,p4) and Orientation(p1,p2,p3) have the same sign
#
m11(p1,p4) = Px[p1]-Px[p4]
m21(p2,p4) = Px[p2]-Px[p4]
m31(p3,p4) = Px[p3]-Px[p4]
m12(p1,p4) = Py[p1]-Py[p4]
m22(p2,p4) = Py[p2]-Py[p4]
m32(p3,p4) = Py[p3]-Py[p4]
m13(p1,p4) = (Px[p1]-Px[p4])**2 + (Py[p1]-Py[p4])**2
m23(p2,p4) = (Px[p2]-Px[p4])**2 + (Py[p2]-Py[p4])**2
m33(p3,p4) = (Px[p3]-Px[p4])**2 + (Py[p3]-Py[p4])**2
Det(p1,p2,p3,p4) = m11(p1,p4)*(m22(p2,p4)*m33(p3,p4) - m32(p3,p4)*m23(p2,p4)) + \
m12(p1,p4)*(m23(p2,p4)*m31(p3,p4) - m33(p3,p4)*m21(p2,p4)) + \
m13(p1,p4)*(m21(p2,p4)*m32(p3,p4) - m31(p3,p4)*m22(p2,p4))
# create triangle data
set print $Triangles
do for [i=1:TriangleCount(0)] {
print sprintf("%g %g",TAx(i,1),TAy(i,1))
print sprintf("%g %g",TAx(i,2),TAy(i,2))
print sprintf("%g %g",TAx(i,3),TAy(i,3))
print sprintf("%g %g",TAx(i,1),TAy(i,1))
print ""
}
unset print
@CreateEdgeDatablock
@CreateHullDatablock
@PlotEverything
set label 1 "Flipping diagonals... "
###
# loop edges and check if need to flip. If on edge was flipped, start over again.
flip = 0
flipCount = 0
flippedAtLeastOnce = 1
while (flippedAtLeastOnce) {
flippedAtLeastOnce=0
do for [e=1:words(EDGES)/2] { # loop all inner edges
# find the 2 triangles with this edge
p1 = EdgeP(e,1)
p2 = EdgeP(e,2)
found = 0
do for [t=1:TriangleCount(0)] { # loop all triangles
tap1 = TAp(t,1)
tap2 = TAp(t,2)
tap3 = TAp(t,3)
p = p1==tap1 ? p2==tap2 ? tap3 : p2==tap3 ? tap2 : 0 : p1==tap2 ? p2==tap3 ? tap1 : 0 : 0
# print sprintf("%d %d %d: %d",tap1,tap2,tap3,p)
if (p!=0) {
if (found==1) {
ta2=t; p4=p;
flip = sgn(Det(p1,p2,p3,p4))==Orientation(p1,p2,p3)
flippedAtLeastOnce = flippedAtLeastOnce || flip
if (flip) {
flipCount = flipCount+1
print sprintf("Flip % 3d: %d-%d with %d-%d",flipCount,p1,p2,p3,p4)
EDGES = changeEdge(e,p3,p4)
TRIANGLES = changeTA(ta1,p1,p3,p4)
TRIANGLES = changeTA(ta2,p2,p3,p4)
@CreateEdgeDatablock
@CreateHullDatablock
@PlotEverything
break # start over again
}
}
if (found==0) { ta1=t; p3=p; found=1}
}
}
}
}
###
# create quadrangles datablock
Center2x(p1,p2) = (Px[p1]+Px[p2])/2. # x-center of 2 points
Center2y(p1,p2) = (Py[p1]+Py[p2])/2. # y-center of 2 points
Center3x(p1,p2,p3) = (Px[p1]+Px[p2]+Px[p3])/3. # x-center between 3 points
Center3y(p1,p2,p3) = (Py[p1]+Py[p2]+Py[p3])/3. # x-center between 3 points
set print $QUADRANGLES
do for [i=1:TriangleCount(0)] {
do for [p=0:2] {
z = Pz[TAp(i,p+1)]
tap1 = TAp(i,p+1)
tap2 = TAp(i,(p+1)%3+1)
tap3 = TAp(i,(p+2)%3+1)
print sprintf("%g %g %g", Px[tap1], Py[tap1], z)
print sprintf("%g %g %g", Center2x(tap1,tap2), Center2y(tap1,tap2), z)
print sprintf("%g %g %g", Center3x(tap1,tap2,tap3), Center3y(tap1,tap2,tap3), z)
print sprintf("%g %g %g", Center2x(tap1,tap3), Center2y(tap1,tap3), z)
print sprintf("%g %g %g", Px[tap1], Py[tap1], z)
print ''
}
}
set print
set label 1 "Coloring areas."
plot $QUADRANGLES u 1:2:3 w filledcurves closed lc palette, \
$EDGES u 1:2:3:4 w vec lw 1.0 lc "grey" nohead, \
$HULL u 1:2:3:4 w vec lw 1.5 lc "blue" nohead, \
$Data u 1:2:3 w p pt 7 ps 0.5 lc palette
### end of code```
Let's condense this down to the nub of the matter:
$QUADRANGLES << EOD
-24.4571 1.79031 -43.7857
-21.0939 1.52621 -43.7857
-17.9859 2.84816 -43.7857
-18.1135 3.64119 -43.7857
-24.4571 1.79031 -43.7857
-17.7308 1.2621 -22.378
-14.7503 3.37709 -22.378
-17.9859 2.84816 -22.378
-21.0939 1.52621 -22.378
-17.7308 1.2621 -22.378
-11.7699 5.49208 -64.641
-18.1135 3.64119 -64.641
-17.9859 2.84816 -64.641
-14.7503 3.37709 -64.641
-11.7699 5.49208 -64.641
-24.4571 1.79031 -43.7857
-21.0939 1.52621 -43.7857
-18.5289 -1.54692 -43.7857
-18.9279 -2.95143 -43.7857
-24.4571 1.79031 -43.7857
-17.7308 1.2621 -22.378
-15.5648 -3.21554 -22.378
-18.5289 -1.54692 -22.378
-21.0939 1.52621 -22.378
-17.7308 1.2621 -22.378
-13.3988 -7.69317 103.08
-18.9279 -2.95143 103.08
-18.5289 -1.54692 103.08
-15.5648 -3.21554 103.08
-13.3988 -7.69317 103.08
-24.4571 1.79031 -43.7857
-19.7591 11.5522 -43.7857
-17.096 9.53213 -43.7857
-18.1135 3.64119 -43.7857
-24.4571 1.79031 -43.7857
-15.0611 21.314 -321.012
-13.4155 13.403 -321.012
-17.096 9.53213 -321.012
-19.7591 11.5522 -321.012
-15.0611 21.314 -321.012
EOD
plot $QUADRANGLES u 1:2:3 with filledcurves closed lc palette
The plot has a border around each quadrangle because the program's default fill style is equivalent to set style fill empty border. The "empty" part is superseded by plot with filledcurves because it makes zero sense to not fill a filled curve. But the "border" part is honored. You can remove the border by changing the global default fill style:
set style fill solid noborder
or you can put this in the plot command instead
plot $QUADRANGLES u 1:2:3 with filledcurves fillstyle solid noborder fc palette
In new-ish gnuplot versions, with filledcurves closed can be replace by with polygons.
For completeness I'll add that even if the fill style does specify a border, you can make the border invisible by using the "nodraw" line type.