nowadays the best way is probably to read the depth position in the vertex shader, depending on the occlusion result, you get transform the vertex properly or set some default value, e.g. (0,0) as position and all flares gonna be culled.

that's minimal latency (occlusion queries are probably just available at least one frame later), minimal cost (it's not some per pixel cost, but just few e.g. 4 vertex read per flare) and low implementation time (just set the depth sampler, z-texture and the proper 2d read position to the lens flare shader).

I read one way, where the effect was drawn whether or not the results of a single pixel depth read back was visible or not (ugh, depth read)

What's wrong with doing depth reads? The lens flares don't need to do depth-testing or depth-writes, so there's no conflict there. Which API are you using?

For the soft variation, I read a 16x16 pixel area around the centre of the light-source, testing for depth values that are closer than the light source.

For each light-source, I draw a 16x16 pixel quad to an intermediate render-target, and the "depth test" shader outputs black or white, depending on whether that pixel is occluded or not.

I then down-sample this intermediate texture 16x, so that I end up with one pixel per light-source, with a fractional occlusion factor (16x16 booleans averaged == 256 possible values).

If I know there's 16 lights that need a fractional occlusion value in the scene, then I'll make my intermediate render target 256x16 pixels (i.e. that's 16 squares, of 16x16 pixels).

The lens flare shader can then sample one value from the downsampled intermediate texture, and get a fractional occlusion value from 0.0f to 1.0f (in 1/256 increments).

Here's a visual example with the scene above, and the intermediate texture below. On the right, the sun's lens flare is still drawn, but with 50% of the regular alpha.

For smaller light-sources (like the individual lights on the stadium spotlight arrays in that picture), I use the method you mention, where the lens-flare itself samples the depth-buffer once, and hides itself if this manual/custom depth-test fails.

You might be able to find a tutorial by searching for "opengl use depth buffer as texture" etc. Being able to use the depth-buffer as a texture is an important ability for shadow-mapping, deferred rendering, etc, so there should be some guides around.
My OpenGL is a bit rusty, but IIRC the process of reading the depth buffer in a shader requires the use of FBO's... You bind a depth-texture to an FBO as the depth/stencil target, you draw your scene, and then afterwards, you can use that depth-texture as if it were a regular old texture. You can then sample from it in a fragment shader, just like you do with any other texture.

the example I was referring to was a single pixel depth read on the cpu (glReadPixels).

Yeah, ok that would be bad

How are you determining the center of the 16x16 area? Manually projecting the light point?
Without using the 16x16 blend method mentioned, or single depth read on cpu compared against the manually projected depth, how would you go about killing
the multiple flare quads based on the depth test result of the main flare point?

Yeah, on the CPU side, you project the light source into screen coordinates. If done correctly, this gives you the x/y position, but also the z value (all these values from -1 to 1 in range). To convert the x/y values to texture coordinates, you can rescale them to the 0-1 range (and maybe flip y upside down if GL's texture coordinates require that).

In the single depth sample method, in your pixel shader, you can sample the depth texture at this x/y position, and discard (or output alpha of zero) if the returned depth value is less than the z value calculated on the CPU.

In the 16x16 method, in your first (intermediate) pixel shader, you sample the depth texture at this x/y position, plus an offset depending on which pixel is being drawn. Then if the returned depth value is less than the z value calculated on the CPU you output black, otherwise output white. Then you calculate mip-maps for this intermediate texture.
Then in the second (actual lens flare) pixel shader, you sample the lowest-resolution mip-map of the intermediate texture, and use the returned value to scale your alpha value (or discard if the sampled value is zero).

If you're using a modernish GPU (e.g. GL3+) then you can do these texture fetches in the vertex-shader instead of for every single pixel in the lens flare.