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.
Description
surface_direction(..., method='geodesic')builds its compass bearing frompixel index offsets scaled by the raster's cell size:
In geodesic mode the coordinates are lat/lon, so
cellsize_xandcellsize_yare 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 theeast-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_xandcellsize_yare already inthe same map unit.
The rest of geodesic mode does the right thing.
_precompute_dd_gridcomputesreal great-circle horizontal distances per pixel and
_dijkstra_geodesicconsumes them, so
surface_distance()andsurface_allocation()are correctin 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 tocompute 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-accuracyreview of #3719. No reproduction scriptattached: 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_directionreports.