I'm making an strategy game using openGL on Android.

It is time to upgrade my unit selection code, which now consists of simple ray to plane intersection.

What I did the last time is use my own math functions and transformations to figure out which unit the user is clicking.

So the question arises: Once I've rendered a scene, is there an easy way to figure out which the unit the user is tapping?

I'm not using any engine here, just openGL on android with my home-brew math library.

The first, revert the clicked screen coordinate to world coordinates by using your projection in reverse. This is done by inverting the modelviewprojection matrix, which if you aren't familiar with it is going to be tricky. Then use the world coordinate to look up the nearest unit.

Why invert?

What I do sounds very similar to what he is already doing: computing the camera matrix (inverse view matrix, without inverting,) this produces the world space position as well as the basis vectors for "camera." When I have window coordinates to pick, I convert them to normalized device coordinates then use the world-space information + projection information (i.e. near plane distance) to construct a ray, then perform a search to find the hit. If scaling is involved in the camera matrix calculation, I use the properties of matrices to avoid normalizing the basis / direction vectors.

I like this approach because I need the world-space camera information for other functions such as culling, billboarding, etc.

If you need pixel perfect picking then the selection buffer approach described by Zoomulator is a simple solution. Though if the user does not tap exactly on a pixel covered by the unit it will not register a hit.

Does anyone have a link to an article, or some sample code?
Doesn't matter which language.

Does anyone have a link to an article, or some sample code?
Doesn't matter which language.

While my personal code is kind of tailored to my application, MathGeoLib's frustum(.h|.cpp) code is very clean (though I personally prefer the "radar" approach to intersection testing.)

An example of what I described:
If your view transform is something like this:
Translate1*Scale*RotateZ*RotateX*RotateY*Translate2

Your camera transform is this:
Translate2' * RotateY' * RotateX' * RotateZ' * Scale' * Translate1'

Now you can extract the position of the camera from the resulting matrix, like in the image on this page.
Note: You'll likely want the left / up / forward vectors to be normalized, what you can do as an optimization is the following:
Translate2' * RotateY' * RotateX' * RotateZ' * Translate1' * Scale' (only apply inverse scaling to the transformation part of the matrix).
That way the columns in the rotation part of the matrix are all unit.

This gives you the camera information.

Take the coordinates of the tapping and convert them from window coordinates to normalized device coordinates, then your ray's direction is simply:
Forward * nearplane distance + Left * x (in NDC) + Up * y (in NDC)
Where forward/left/up are extracted from the camera matrix as described above, the ray starts from your camera's position (translation part extracted from the matrix.)
Note, you might have to change some of the signs for the direction equation depending on the convention you use for the window coordinates origin.

Disclaimer: I don't remember why there's a conversion to NDC nor do I have the time to verify it right now.
Also, I had written in my notes, that this assumes a symmetric frustum, but I can't remember why either.

There's another alternative, less flexible but should work nonetheless. Make a 3d cursor like you would make any other unit, as a 3d model hovering the terrain, and its coordinates on the terrain are controlled by the mouse. It won't work with the UI, though...

This runs on a touch screen. Hence, there is no mousemove event before the screen is tapped.