To have "map scale factors [...] included implicitly," taking inspiration from #1389 (comment), I think this could be implemented something like

def jacobian(field_0, field_1):
    grad_field_1_y, grad_field_1_x = gradient(field_1, axes=('y', 'x'))
    return vorticity(field_0 * grad_field_1_x, field_0 * grad_field_1_y)

Since

Our current implementation of vorticity won't work with @jthielen suggestion since we check for units of speed in a decorator with the function. Doing two gradient calls doesn't seem like it would be too much with a function like

def jacobian(field_0, field_1):
    grad_field_1_y, grad_field_1_x = gradient(field_1, axes=('y', 'x'))
    grad_field_0_y, grad_field_0_x = gradient(field_0, axes=('y', 'x'))
    return grad_field_0_x * grad_field_1_y - grad_field_0_y * grad_field_1_x

Otherwise, we could develop a more generic curl function that might help both vorticity and this calculation.

@kgoebber That approach was my first thought, but I avoided it because my prior impression from #893 was that scalar gradients (as used there) wouldn't have any scale factor corrections, and so this in turn would lack them, compared to the curl-based approach which would.

However, this mismatch between formulations that should otherwise be equivalent worried me a bit, so I went back to https://en.wikipedia.org/wiki/Orthogonal_coordinates and the implementations in #893 to figure out the proper/expected behavior. Long story short, I'm pretty sure we need to apply scale factor corrections in gradient:

so that the Jacobian (via either curl of product with gradient or cross product of gradients, the former results in a lot of cancellations so the latter is probably better computationally) ends up as:

Ah, yes. Good catch.