Skip to content

Saved coordinate values are different with numpy2. #686

@pp-mo

Description

@pp-mo

Coordinate values changed in a saved grib file, when numpy v1.26.4 was replaced by v2.2.6

The example has coordinates are derived from a PP file, with the following original header word values:

FLD:  lbrow 1920
FLD:  lbnpt 2560
FLD: bzx -0.070312500000000000000000000000
FLD: bdx 0.140625000000000000000000000000
FLD: bzy -90.046875000000000000000000000000
FLD: bdy 0.093750000000000000000000000000

this loads as XY coordinates (unbounded) with the following endpoints:

Lons: 0.0703125 ... 359.9297
Lats: -89.953125 ... 89.953125

and the following step values:

Lon step value: 0.140625
Lat step value: 0.09375

The save code converts these to these key values in the saved file:

  • iDirectionIncrement = 140625
  • jDirectionIncrement = 93750
  • longitudeOfFirstGridPoint = 70312
  • longitudeOfLastGridPoint = 359929687 (numpy 1.26) --> 359929696 (numpy 2.2)
  • latitudeOfFirstGridPoint = -89953125 (numpy 1.26) --> -89953128 (numpy 2.2)
  • latitudeOfLastGridPoint = 89953125  (numpy 1.26) --> 89953128 (numpy 2.2)

Cause + simple fix

The reason for the differences is that numpy2 doesn't promote value types in the same way.
The actual calculation is here
Notably, this does effectively int(value * 1000000).
But, numpy2 reproduces the previous values if we instead use int(value.astype("f8") * 1000000)

I think that the latter is equivalent to what happened in numpy 1, i.e. when multiplying by an integer the float32 was previously automatically promoted to a float64 (?maybe just by conversion to Python float and back again?). So making that explicit delivers the same result using numpy 2.

Additional Requirement

However the original different values caused an error in other software, which was apparently checking that, e.g.
longitudeOfLastGridPoint == longitudeOfFirstGridPoint + iDirectionIncrement * (Ni - 1)
The above calculation fix doesn't necessarily guarantee that, especially as our source is 32-bit PP data, which has caused us previous problems with the rounding of coordinate values.

So possibly, the fix should ensure the expected relationship between (start, stop, step, number), when these are derived (by scaling + rounding) from the coordinate point values.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    Status

    No status

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions