Why does ReliefPlot not have a version for functions?

TL, DR: Use the LightingAngle option to DensityPlot to achieve the same shadows effect as ReliefPlot. (@Brett Champion's answer)

To plot bivariate functions, we typically either use Plot3D, whose dynamic effects allow flexible investigation, or DensityPlot, which is better for publications. However, sometimes the output of DensityPlot is a little ... obscure(?), and shadows from ReliefPlot would add a lot perspective.

Example from the ReliefPlot documentation:

ReliefPlot[Table[i + Sin[i^2 + j^2], {i, -4, 4, .03}, {j, -4, 4, .03}], ColorFunction -> "SunsetColors"]


Compare a replication using DensityPlot:

DensityPlot[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4}, ColorFunction -> "SunsetColors", PlotPoints -> 40]


Where the well-known oddity of ReliefPlot's orientation can be seen.

For years on MMA SE, a lot of members have used ReliefPlot via manual sampling for function plotting questions (e.g., this post), yielding better aesthetic results.
Alternatives include:

• Using Plot3D and set the viewpoint at infinity above, but it can also hairy to do. (See @MichaelE2's answer which does this elegantly)
• Manually define normal vectors and do a normal mapping, which I am not sure how to implement in Mathematica right now. (See my answer to get a crude idea)

But why isn't ReliefPlot for functions been officially implemented (i.e., with adaptive sampling like other *Plot* built-ins)? What is the real technical difficulty?

P.S.: I believe extensive about this problem can already be found on MMA SE (like here and here), but they are scattered everywhere so I am raising this question as a port.

For this thread to serve as a guide to others with the same problem, what are the ways to mimick ReliefPlot's texture? (please provide simple code and output examples) From your personal experience, at what point should one bother to do this for the nicer touch?

You can use the LightingAngle option to DensityPlot to achieve the same shadows effect as ReliefPlot:

DensityPlot[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4},
ColorFunction -> "SunsetColors", PlotPoints -> 40,
LightingAngle -> 180*Degree]


• Interesting, didn't know this option. According to the document of LightingAngle, the default setting of ReliefPlot seems to be LightingAngle -> 3/4 Pi? Commented May 16, 2021 at 3:05
• Yes, that's correct. For my example I just used the first suggested completion for the option. Commented May 16, 2021 at 14:42
• Neato! I wonder if this option was available all along in recent versions of MMA... Looked back at the documentation, wondering why I missed this, and find that LightingAngle is listed in the header section, but is lacking from the examples section. Commented May 18, 2021 at 8:52
• @Gravifer It's already there since v7: i.sstatic.net/e4Ov4.png Strangely no example is added to the document of DenistyPlot. Commented May 19, 2021 at 3:08

ReliefPlot is using the 3D lighting model (I think) to highlight changes in the gradient.

Deploy@Plot3D[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4},
ColorFunction -> (ColorData["SunsetColors"][#3] &),
PlotPoints -> 35, MaxRecursion -> 4, Mesh -> None,
ViewPoint -> {0, 0, Infinity},
Lighting -> {{"Ambient", GrayLevel[0.75]}, {"Directional",
GrayLevel[0.15], ImageScaled[{-2, 0, 2}]}, {"Directional",
GrayLevel[0.15], ImageScaled[{-2, -2, 2}]}, {"Directional",
GrayLevel[0.15], ImageScaled[{0, -2, 2}]}}, Boxed -> False,
Axes -> False]

• Really elegant results. But I do wonder if ReliefPlot actually constructs a 3D surface internally; it is a function often used on large geographic datasets, so actual ray tracing would be quite expensive. I do wonder that if the plot is complicated (e.g., when poles occur), employing Plot3D may lead to performance issues. I suspect ReliefPlot calculates the gradient of the array and do a normal mapping accordingly. Commented Mar 9, 2021 at 4:52
• @Gravifer Well, yeah, of course. Plot3D does not do raytracing either, but calculates the normals (symbolically or numerically, depending on the function). Graphics3D uses the normals (not Plot3D) and Lighting to determine the colors of the polygons. (Unless that is done by the GPU. I'm not really up on exactly what GPUs are used to do.) ReliefPlot constructs a Raster of the colors that it determines. Commented Mar 9, 2021 at 4:55
• I am a little confused. I just verified that ReliefPlot returns a Graphics object; does it calls upon Graphics3D functionalities internally? I didn't realize that it uses polygon based reflection rendering. Commented Mar 9, 2021 at 5:12
• I don't know how ReliefPlot actually does what it does internally, but only inferred it from how well it matches the Plot3D output. It probably uses the same basic algorithm for 3D shading to determine how to color each square in the raster, but it might use a routine optimized for producing a 2D raster, instead of rasterizing a 3D plot. Note that ReliefPlot has an option LightingAngle. That and other language in the docs suggest that it produces a 2D rendering of a 3D surface. The Method option can be used to switching between two shading algorithms. Commented Mar 9, 2021 at 5:50
• @Gravifer This might give a slightly closer reproduction, though the shadow transitions are still not a diffuse: Plot3D[i + Sin[i^2 + j^2], {i, -4, 4}, {j, -4, 4}, ColorFunction -> (ColorData["SunsetColors"][#3] &), PlotPoints -> 265, MaxRecursion -> 0, Mesh -> None, ViewPoint -> {0, 0, Infinity}, Lighting -> {{"Ambient", GrayLevel[0.85]}, {"Directional", GrayLevel[0.25], ImageScaled[{-1, 1, Sqrt[2]}]}}, Boxed -> False, Axes -> False] Commented Mar 9, 2021 at 14:42

Update Found this ResourceFunction; really nice touch.

As this is a guideline question, I'll place my own research in a standalone answer.
The ideas comes from the discussion with @MichaelE2 under his answer following the Plot3D approach.

I mentioned there that it is possible do do artificial normal mapping and adding the shadow as a mask over the original DensityPlot. The following code is an example for how to do this:

theta=3Pi/4; phi=Pi/4;
func = j + Sin[i^2 + j^2];
base = DensityPlot[func, {i,-4, 4}, {j,-4, 4}, ColorFunction->"SunsetColors", PlotPoints->40];
DensityPlot[
Evaluate[
FromSphericalCoordinates[{1,theta,Pi-phi}](*Default source position of ReliefPlot*)
. Normalize@Append[-Grad[func, {i, j}], 10](*2nd argument of GrayLevel is Opacity*)
], {i,-4, 4}, {j,-4, 4}, ColorFunction->(GrayLevel[#,(2#-1)^2]&), PlotPoints->80],
Opacity[.5], {1, 1}];
ReliefPlot[Table[i + Sin[i^2 + j^2],
{i,-4, 4, .03}, {j,-4, 4, .03}], ColorFunction->"SunsetColors"]}]


Note that:
(1) I switched i and j so the orientation matches.
(2) The $$z$$ component I appended to the -Grad vector is $$10$$ instead of $$1$$; this effectively shrinked the surface and smoothed the shadows. You may tweak this to achieve better visual effects.

The output is still isn't really quite the same to the built-in, but can show the gist of it. I used a quadratic opacity scaling in the mask; it can be modified to improve the output. What you need to know is that in a true normal mapping, the shadow mask is multiplicative; a additive mask will never create highlight faithfully.

Note that there are a few artifacts: the direction of the normal is somewhat 'singular' where the gradient vanishes.

Pack this to a single function (can be really buggy)

reliefPlot[f_, {x_, xmin_, xmax_}, {y_, ymin_, ymax_},
ops : OptionsPattern[{LightingAngle -> {(3 \[Pi])/4, \[Pi]/4},
DensityPlot, ReliefPlot}]] :=
Show[{DensityPlot[f, {x, xmin, xmax}, {y, ymin, ymax}, ops],
Insert[DensityPlot[
Evaluate[
FromSphericalCoordinates[{1,
OptionValue[LightingAngle][[1]], \[Pi] -
OptionValue[LightingAngle][[-1]]}] .
Normalize@Append[-Grad[f, {x, y}], 10]], {x, xmin, xmax}, {y,
ymin, ymax}, ColorFunction -> (GrayLevel[#, (2 # - 1)^2] &),
Evaluate[Sequence @@ FilterRules[{ops}, PlotPoints]]],
Opacity[.5], {1, 1}]}]