I’m creating a 3D globe with elevation:

I now want to add click detection to get a 3D world xyz position from a 2D screen xy + depth (from a click). I’m using the following code (inspired from Getting World Position from Depth Buffer Value) to unproject the 2D coordinate + depth (using gl-matrix for computations) into 3D space:

static unProject(screenX: number, screenY: number, depth: number, projMatrix: mat4, viewMatrix: mat4): vec3 {

  // Convert x, y, depth from [0 .. 1] to NDC [-1 .. 1]
  const clipSpacePosition: vec4 = vec4.fromValues(
    screenX * 2.0 - 1.0,
    screenY * 2.0 - 1.0,
    depth * 2.0 - 1.0,
    1.0
  );

  let viewSpacePosition: vec4 = vec4.transformMat4(vec4.create(), clipSpacePosition, mat4.invert(mat4.create(), projMatrix));

  // Normalize eye coordinates
  viewSpacePosition[0] /= viewSpacePosition[3];
  viewSpacePosition[1] /= viewSpacePosition[3];
  viewSpacePosition[2] /= viewSpacePosition[3];
  viewSpacePosition[3] /= viewSpacePosition[3];

  const worldSpacePosition: vec4 = vec4.transformMat4(vec4.create(), viewSpacePosition, mat4.invert(mat4.create(), viewMatrix));

  // Normalize world coordinates
  worldSpacePosition[0] /= worldSpacePosition[3];
  worldSpacePosition[1] /= worldSpacePosition[3];
  worldSpacePosition[2] /= worldSpacePosition[3];

  return vec3.fromValues(worldSpacePosition[0], worldSpacePosition[1], worldSpacePosition[2]);
}

Using this function I clicked on the south-east of France and plotted 1000 cubes having a depth in the [0 .. 1] range:

Without the sphere:

Cubes seems correctly displayed between my z-near & z-far planes, but the resolution near z-near is much higher than near z-far (cubes are not evenly distributed along the axis). If I reduce z-near, it worsen this effect (and at the extreme all the cubes are drawn on z-near). This is normal according to https://learnopengl.com/Advanced-OpenGL/Depth-testing , but may indicate a problem?

Empirically I found the depth value of the cube intersecting my sphere (the one in purple). I now need to get it from OpenGL depth buffer. I thus drawn my scene again with a simplified fragment shader encoding the depth (gl_FragCoord.z) into RGBA:

#version 300 es
precision highp float;

vec4 encode(float depth) {
    vec4 bitSh = vec4(256.0 * 256.0 * 256.0, 256.0 * 256.0, 256.0, 1.0);
    vec4 bitMsk = vec4(0.0, 1.0 / 256.0, 1.0 / 256.0, 1.0 / 256.0);
    vec4 enc = fract(depth * bitSh);
    enc -= enc.xxyz * bitMsk;
    return enc;
}

void main(void) {
    color = encode(gl_FragCoord.z);
}

And decoded the value on JS side this way:

decode(color: vec4) {
  const bitSh: vec4 = vec4.fromValues(1.0 / ( 256.0 * 256.0 * 256.0 ), 1.0 / ( 256.0 * 256.0 ), 1.0 / 256.0, 1.0);
  return vec4.dot(color, bitSh);
}

I also written the encode method in JS to ensure that encoding / decoding a value is working, it does (with a little precision error, for instance if I encode / decode 0.962 I get 0.9619999999999997).

Using this method, the depth value retrieved from OpenGL is too small (in my example with a z-near of 1 & z-far of 500, I’m expecting ~0.962 and I have ~0.268).

Am I missing a transformation on my depth buffer value?