几天后，我找到了适合我的解决方案。也许有人会发现它有用：

module PrintToGraph where

import qualified Graphics.Gnuplot.Advanced as GP
import qualified Graphics.Gnuplot.Frame as Frame
import qualified Graphics.Gnuplot.Frame.OptionSet as OptsSet
import qualified Graphics.Gnuplot.Plot.ThreeDimensional as Plot3D
import qualified Graphics.Gnuplot.Graph.ThreeDimensional as Graph3D
import qualified Graphics.Gnuplot.LineSpecification as LineSpec
import GHC.Exts (groupWith )
import qualified Graphics.Gnuplot.Value.Atom as Atom
import Graphics.Gnuplot.ColorSpecification ( paletteFrac )
import Data.Foldable ( Foldable(foldMap') )
import Data.List ( elemIndex )
import Data.Maybe ( fromJust )

    defltOpts :: OptsSet.T (Graph3D.T Double Double Double)
defltOpts = OptsSet.key False  OptsSet.deflt

waveFuncVis :: (Double -> (Double, Double, Double) ->  Double) -> Double -> Double -> Frame.T (Graph3D.T Double Double Double)
waveFuncVis func depth precision =
      let x = Plot3D.linearScale 100 (-10, 10)
          testedRange = (groupWith (\(x,y,z) -> test func (x,y,z) depth precision)  . filter (\(x,y,z) -> funcWrapper func x y z^2 >= precision)) [(x1,y1,z1) | x1<-x, y1<-x, z1<-x]
          range = [(x1,y1,z1) | x1<-x, y1<-x, z1<-x]
          calcColor :: [(Double,Double,Double)] -> Double
          calcColor array  =  fromIntegral (fromJust (elemIndex array testedRange)) / fromIntegral (length testedRange)
          linespec array   = Graph3D.lineSpec $ LineSpec.lineColor (paletteFrac (calcColor array)) LineSpec.deflt
          graph array  = linespec array <$> Plot3D.cloud Graph3D.points array
    in  Frame.cons defltOpts $ foldMap' graph testedRange

test :: (Double -> (Double, Double, Double) -> Double)
    -> (Double, Double, Double) -> Double -> Double -> Integer
test func (x, y , z) depth precision
    |  funcWrapper func x y z^2 >= precision = round $ funcWrapper func x y z^2 * depth
    |  otherwise = 0

funcWrapper :: (Double -> (Double, Double, Double) ->  Double) -> Double -> Double -> Double -> Double
funcWrapper func x' y' z' = func 1.0 (toR x' y' z', toTau x' y' z', toPhi x' y' z')

--2pz Hydrogen function
waveHfunc2pz :: Double -> (Double, Double, Double) ->  Double
waveHfunc2pz z (r, tau, phi) = a * b * c* e
                            where a,b,c,e :: Double
                                  a = 1.0/(4.0*sqrt (2.0*pi))
                                  b = (z/aBohr)**2.5
                                  c = pureTrig cos tau
                                  e = r*exp(-1.0 * (z*r/(2.0*aBohr)))

main ::  IO ()
main = sequence_ [GP.plotDefault (waveFuncVis waveHfunc2pz 10000 0.0005)]

简要地说：
1.我们丢弃小于 precision 的函数值。（为此，我在testedRange中使用了filter）
1.由于groupWith，我们得到了坐标列表的列表-（x，y，z）。这里的每个子列表都包含了给出相同函数值的坐标。
1.为了给它们着色，我们将子列表的索引转换为Double值，并将其用作PaletteFrac的参数。
结果，我们得到一团彩色点，其中每种颜色对应于一个函数值。
Example图像。