MetOffice / MetOffice/CSET

Data is plotted offset to coastlines due to PROJ 9.8

Open
#2,052 5 comments 1 reaction 0 assignees View on GitHub
bug
Dominant language
Python
Stars
33
Forks
19
Avg merge
1d 22h
Merged PRs (30d)
30

Description

### Describe the bug

When using PROJ 9.8 the data is offset from the coastlines/axis labels when plotting, as can be seen here:

Topography data for the South West UK, with the data shifted to the south of where it should be.

The issue was occurring when plotting from within the CSET environment, while not from older environments, such as the scientific software stack. @ukmo-huw-lewis produced a test script to easily demonstrate the issue, which proved invaluable in narrowing it down. Given it was environment dependent, and could be reproduced without any CSET code running, it was definitely an environment issue.

```py
import iris
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

fname = "/path/omitted/.../umnsaa_pvera000"
cube = iris.load_cube(fname, "surface_temperature")[0]
coord_system = cube.coord("grid_latitude").coord_system
projection = ccrs.RotatedPole(
pole_longitude=coord_system.grid_north_pole_longitude,
pole_latitude=coord_system.grid_north_pole_latitude,
central_rotated_longitude=0.0,
)
figure = plt.figure()
ax = figure.add_subplot(projection=projection)
clat = cube.coord("grid_latitude").points
clon = cube.coord("grid_longitude").points
ax.pcolormesh(clon, clat, cube.data)
ax.coastlines()
ax.gridlines(alpha=0.3, draw_labels=True)
plt.show()

```

I took this chance to finally learn [git bisect](https://git-scm.com/docs/git-bisect), which allows for efficiently searching through commits to find when something changed/broke.

The result of my bisecting identified commit 0dd92466e1a1a310b0101b1efaa17a0009910cb1 (9th of March) as the one that caused the breakage. Unfortunately it was quite a big environment update, so identifying the specific package took a little more work, eventually turning out to be PROJ.

The end result of this investigation is that the issue was coming from the PROJ 9.8.0 release, which included some geodesic changes that change results for some common transforms. I’ll be honest, the details go over my head a little, but there is a good discussion of the issues in https://github.com/OSGeo/PROJ/pull/4736.

The 9.8.0 package on conda-forge has been marked as broken so when the lockfile update action runs it will either go back to 9.7.1, or if 9.8.1 has been released by then it will use that. By the time Monday rolled around [PROJ v9.8.1](https://github.com/OSGeo/PROJ/releases/tag/9.8.1) had been released, and that was what fixed this issue in https://github.com/MetOffice/CSET/pull/2043#pullrequestreview-4097351542.

Contributor guide

No contributing guide indexed for this repository

Research direction

Start with the reproduction script in the issue and inspect the linked pull request #2043 to follow the dependency change that addressed the offset. Verify the plotting example with the relevant environment and confirm that the data aligns with coastlines and axis labels.

Written by the indexing model from the issue text.

Assessment

Tech stack
python
Domain
data-visualization
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Stale
Clarity
Clearly specified
Newbie friendliness
15/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.