How to create and export a procedural planet texture with normal map

Here is something to get you started. Go to Help > Compile > Neat Examples to get code for Perlin noise. Make an image with, for example,

amplitude = 1.5;
frequency = 1.0;
gain = 0.5;
lacunarity = 2.0;
increment = 7.0;
scale = 150.0;
sphereImage = ReliefImage[perlin[512, 512]]

Several answers here have described how to map an image to a sphere. For example, here and here. Unfortunately, artefacts seem unavoidable; perhaps someone can suggest a better way.

ParametricPlot3D[
   {Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]},
   {u, 0, 2 Pi}, {v, 0, Pi},
   Mesh -> None, PlotPoints -> 100, TextureCoordinateFunction -> ({#4, 1 - #5} &),
   Boxed -> False, PlotStyle -> Texture[Show[sphereImage, ImageSize -> 400]], 
   Lighting -> "Neutral", Axes -> False]

SphericalPlot3D[1, {theta, 0, Pi}, {phi, 0, 2 Pi},
   Mesh -> None, 
   TextureCoordinateFunction -> ({1 - #4, #5} &), 
   PlotStyle -> Directive[Texture[sphereImage]], Lighting -> "Neutral", 
   Axes -> False, Boxed -> False]

POVRay uses 3D Perlin noise, and generates much more realistic planets than my examples here. Plus, there are no artefacts.

One way to approach the generation of a planetary surface is to start with something that is "like" what you want and then extrapolate beyond the boundaries. With images, this is sometimes called "inpainting" and Mathematica has the function Inpaint (which, according to the documentation, was introduced in version 8). For instance, here is an imaginary Earth (downloaded from google images). A small portion of this is selected and then a larger map is extraplated from the selection:

img = Import["http://www.geocurrents.info/wp-content/uploads/2012/01/Map_of_\
Earth_2.jpg"]; 
imgSmall = ImageCrop[img, {150, 150}]
dims = ImageDimensions[imgSmall];
canvas = ImageCrop[imgSmall, 2*dims, Padding -> White];
mask = ImageCrop[ConstantImage[Black, dims], 2*dims, Padding -> White];
extrapolated = Inpaint[canvas, mask, 
  Method -> {"TextureSynthesis", "NeighborCount" -> 20, "MaxSamples" -> 1000}]

And here is the same using the picture of Mars in the OPs link:

For mapping to the sphere, you can:

colorfun = BSplineFunction[ImageData[extrapolated], SplineDegree -> 1];
ParametricPlot3D[{Sin[u] Cos[v], Cos[u] Cos[v], Sin[v]}, {u, 0, Pi}, {v, -Pi, Pi}, 
  ColorFunction -> {colorfun[1 - #5, #4] &}, ViewPoint -> {Pi/2, Pi, 1}, 
  PlotPoints -> 200, MaxRecursion -> 0, Mesh -> None, Axes -> None, Boxed -> False, 
  BoundaryStyle -> None, ImageSize -> {500}]

Being a 3D plot, you can rotate the planet at will. There is an unfortunate stitching artifact using this mapping, but for this image it's on the other side of the planet.

I'm fairly certain that the OP has move on to other projects. That said, the other answers don't address one part of the question, which is the normal map generation for inclusion in 3D software.

Let's take an example texture and use the following function to generate a normal map.

imgd = ExampleData[{"ColorTexture", "WhiteMarble"}];
makeNormalMap[i_Image] := 
 Module[{grayscale, norm, didx, didy, nx, ny, nz},
  grayscale = ColorConvert[i, "Grayscale"];
  didx = ImageData@DerivativeFilter[grayscale, {0, 1}];
  didy = ImageData@DerivativeFilter[grayscale, {1, 0}];
  norm = Sqrt[didx^2 + didy^2 + 1];
  nx = Rescale[-didx/norm];
  ny = Rescale[-didy/norm];
  nz = 1/norm;
  Image[{nx, ny, nz}, Interleaving -> False]
]

The process for making the normal map is adapted from a similar question at Blender Stackexchange. Note: I would be interested in some code-review here, since I suspect my hackish implementation of Norming ImageData can be improved.

The two images from Mathematica:

For a 3D material example, I am using Unreal Engine to create a material, applying the left-hand image as the base-color texture and the right-hand image as the normal map. Assigning arbitrary specular and roughness parameters (very specular and very smooth), I get the following:

And just for comparison, here is the same material, but without using the normal map.