This issue is ultimately due to #8836. You can follow the arguments there as to why the change was made. As far as I am concerned I do not believe these changes should be reverted (but would be happy to hear why they should, taking into account the arguments linked above).

A nan in one axis will propagate to the other axis because transData relies on matrix multiplication, and the semantics of IEEE floats are such that

[ a 0 ]
[ 0 b ]

times

[ nan ]
[  c  ]

(where a, b and c are finite) will leak the nan into the y axis (to be fully complete, we have a third column in the matrix for the offset, but it does not matter here). I don't think there's much that can be done here.

To restore previous behavior use ax.set_yscale('log', nonposy='clip').

I see your point that this would a result of matrix multiplication. (Can it actually be used for log axis?) But if axes are separable, use of a matrix multiplication seems unnecessary overkill. Gives you two extra superfluous multiply-add with zero (one per axis). And can cause trouble, as we have seen. For this matrix multiplication a non-IEEE multiplication where 0 * nan = 0 would be preferable ...

In my application I define an x-scale that does not really know anything about the y-scale - and should not have to. Hence the trouble. My scale only wants to convert x-coordinates, and not know anything about y coordinates ... A solution could be if transforms that are separable have methods transform_x and transform_y to only transform those coordinates w/o carrying about the other. Maybe it already exists?

Though this would be even more overkill, the technically correct solution for separable axes might be to use a sparse matrix. Then you never multiply the nan with a zero and add it to anything.

https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.get_xaxis_transform.html#matplotlib.axes.Axes.get_xaxis_transform is what you want. See also https://matplotlib.org/users/transforms_tutorial.html#blended-transformations.
Arguably that section of the transforms tutorial could use some rewording. If you have any suggestions for that (given that you are actually using the feature) feel free to suggest them.

The linear part of the transform is done at C-level by the rendering code, so using anything but a C-style array to represent the array is going to be tricky at best.

Yes, this does the trick. Thank you!

import matplotlib.pylab as plt
f = plt.figure()
ax = f.add_subplot(111)
ax.set_yscale('log')
ax.get_xaxis_transform().transform([0,0])

Output

array([ 81.5 ,   53.79])

One point I have been struggling with is to get the transform from data to axes coordinates. I have been using (ax.transData + ax.transAxes.inverted()) (now (ax.get_xaxis_transform() + ax.transAxes.inverted()) for just the xaxis values) but surely there must be a better way? In the transform tutorial, there is several formula provided for transData, depending on situation, so one cannot just take the parts, e.g.,transScale and transAxes because in case of projections the setup is different.

No, everything is rooted onto the display (pixel) coordinates so you should just go through them. Well perhaps it would have been nice if {transProjection+transProjectionAffine} == transLimits, but that's probably quite low in the priority list of refactors (as you have noted yourself it's quite easy to break things when fiddling with this kind of internals :-))

Closed by #9375, but don't let the conversation stop because of that!