Geodesic surface_direction scales index offsets by degrees, ignoring meridian convergence

Description

surface_direction(..., method='geodesic') builds its compass bearing from
pixel index offsets scaled by the raster's cell size:

# xrspatial/surface_distance.py, _finalize_direction
dx = (src_col.astype(np.float64) - col_idx) * cellsize_x
dy = (src_row.astype(np.float64) - row_idx) * cellsize_y

In geodesic mode the coordinates are lat/lon, so cellsize_x and cellsize_y
are both in degrees. A degree of longitude is shorter than a degree of
latitude everywhere except the equator, by a factor of cos(lat), so the
east-west leg of that offset is overstated and the bearing tilts toward east
or west. At 60 degrees latitude the longitude degree is half the latitude
degree, and a source one cell northeast reads as roughly 63 degrees rather
than 45.

The planar path is fine: there cellsize_x and cellsize_y are already in
the same map unit.

The rest of geodesic mode does the right thing. _precompute_dd_grid computes
real great-circle horizontal distances per pixel and _dijkstra_geodesic
consumes them, so surface_distance() and surface_allocation() are correct
in geodesic mode. Only the bearing is affected.

This is the residual half of #3719, which fixed the sign of the offset but not
its scale. Splitting it out because the fix needs a decision the sign fix did
not: whether to scale the longitude leg by cos(lat) at the pixel, or to
compute a proper forward azimuth from the two lat/lon pairs, which is the
geodesically correct answer but changes what "bearing" means for long paths
(the initial azimuth of a great circle is not the azimuth at its far end).

Noticed during the /sweep-accuracy review of #3719. No reproduction script
attached: geodesic mode is numpy-only and the defect is a scale factor in a
formula rather than a backend divergence, so the code reads as the evidence.

Expected behaviour

Geodesic bearings should account for the convergence of the meridians, and the
docstring should say which azimuth convention surface_direction reports.

Start in xrspatial/surface_distance.py at _finalize_direction, then read _precompute_dd_grid and _dijkstra_geodesic to understand the existing geodesic path. Resolve whether bearing uses local convergence scaling or forward azimuth, and document the chosen azimuth convention. Done means geodesic direction accounts for meridian convergence while planar behavior remains unchanged.