Coordinates are scrambled for areas crossing the zero meridian

[yellowcap]

What happened?

When using feature crops in polytope, there are two problems in the data returned

1. X and Y are swapped.
2. longitude values are above 360

When requesting data from polytope using a feature that crosses the zero meridian, the data is returned with coordinates that have values above 360. It looks like all points that have longitude above 0 get somehow pushed above 360.

What are the steps to reproduce the bug?

This script loads data, then substracts the values that are too high and swaps X and Y, which results in the plot as expected.

import earthkit.data
import geopandas as gpd
import shapely
import matplotlib.pyplot as plt


request = {
    "class": "d1",
    "dataset": "climate-dt",
    "activity": "ScenarioMIP",
    "experiment": "SSP3-7.0",
    "model": "IFS-NEMO",
    "generation": "1",
    "realization": "1",
    "resolution": "high",
    "expver": "0001",
    "stream": "clte",
    "time": "0000",
    "date": "20260101",
    "type": "fc",
    "levtype": "o2d",
    "param": "263101",
    "feature" : {
        "type" : "boundingbox",
        "points" : [[-2, -1], [2, 1]],
    },
}


ds = earthkit.data.from_source("polytope", "destination-earth", request, stream=False, address="polytope.lumi.apps.dte.destination-earth.eu")

da = ds.to_xarray()

df = da.to_dataframe()
# Transform the data so that the coords are fixed
df.loc[df.y > 180, "y"] = df.loc[df.y > 180, "y"] - 360
# X and Y are swapped in the data
geometry = shapely.points(df.pop("y"), y=df.pop("x"))

gdf = gpd.GeoDataFrame(df, geometry=geometry)
gdf.shape
gdf.plot(
    column='avg_tos',
    cmap='viridis',
)
plt.show()

Renders

Version

v0.11.1

Platform (OS and architecture)

osx

Relevant log output

Accompanying data

No response

Organisation

Development Seed

[mathleur]

Hi @yellowcap,

To swap x and y, please use the following syntax in the feature definition, which explicitly defines the order of the axes to return:

"feature": {
        "type": "boundingbox",
        "points": [[-2, -1], [2, 1]],
        "axes": ['longitude', 'latitude']
    },

For values returned in the longitude range [0, 360], this is returning the coordinates of the data points as stored in the datacube. The returned data points seem to be the right ones however, only translated within the range of the stored data.

Since there is this translation, some post-processing is needed to get it into your desired ranges for plotting. This post-processing step can also be done automatically by earthkit-plots if this might be easier/of interest through the following code:

da = ds.to_xarray()
import earthkit.plots
chart = earthkit.plots.Map()
chart.point_cloud(da['avg_tos'], x="x", y="y")
 
chart.coastlines()
chart.borders()
chart.gridlines()
 
chart.title("{variable_name} (number={number})")
 
chart.legend()
 
chart.show()

I hope this helps answer your questions!

[yellowcap]

hi @mathleur thanks for your response.

The axes fix does not help, the coordinates were already passed in lon/lat order. So adding the long/lat axes indication is not changing the data, x and y are still swapped when I add that. Latitude should map to y, and longitude to x, which is not the case. Its still the other way around.

Regarding the conversion of the values, I did notice that earthkit is able to plot the values correctly out of the box. Why not doing that conversion right away so that the values are in a standard range for the users after download? To me this was very confusing and it took me quite a while to figure out how to transform the data so that it is in correct ranges.

[awarde96]

Hi, in regards to the axes. I have tested and on my end

"feature": {
        "type": "boundingbox",
        "points": [[-2, -1], [2, 1]],
        "axes": ['longitude', 'latitude']
    },

and

"feature": {
        "type": "boundingbox",
        "points": [[-2, -1], [2, 1]],
        "axes": ['latitude', 'longitude']
    },

Do swap the x and y axes.

By default, with no axes we assume "axes": ['latitude', 'longitude']

I am also looking at deploying a branch ecmwf/covjsonkit#60 that returns the names of the axes e.g. latitude longitude in the covjson and xarray so hopefully that will resolve the axes issue and you can use those instead of x and y. I just need to update examples and documentation in quite a few areas before deploying

[yellowcap]

Hi, thank you very much four all the feedback. I tried it and for bbox everything seems to work as expected indeed.

In my real use case I used polygon features and not bounding box, and there the axes seem to be swapped. In the example shared with you I simplified this to use bbox, and in that one it works. Could it be a problem of the polygon feature request?

Could you try the following and see if it works for you?

    "feature" : {
        "type" : "polygon",
        "shape" : [
            [-2, -1],
            [-2, 1],
            [2, 1],
            [2, -1],
            [-2, -1]
        ],
        "axes": ['longitude', 'latitude'],
    },

In this case x and y values get mixed up on my system.

Also, any input on the strange values >360? Do you know where that comes from and how earthkit works around it?

[mathleur]

Hi @yellowcap ,

both the polygon and bounding box features seem to swap the axes in my tests.

I'm also unable to reproduce the values > 360 problem. Would it be possible to know when and how you obtain them? All longitude values returned for me are in the range [0, 360] as expected

[yellowcap]

Hi @mathleur , for polygon, I can get the axes to swap using the axes parameter, but the values are the wrong way around on my system. If I run it as is above, the x values should be in [-2, 2], but they are between [-1, 1], so the latitude maps to x instaead of y.

Also for the values, yes apologies again, the values are in range [0, 360], but that is not as expected for me. In my expectations they would be between [-180, 180]? But that might be because I am not a weather specialist. Is it standard in weather models to have longitude to range [0, 360]? If that is normal, then the transformation below makes sense to me and that is resolved.

df.loc[df.y > 180, "y"] = df.loc[df.y > 180, "y"] - 360

[mathleur]

Hi @yellowcap,

at the moment, our standard is that latitude maps to x and longitude to y, but as @awarde96 mentioned, we will soon be deploying a branch that changes the axis names to latitude and longitude.

In our models, the standard for the range is [0, 360] and I think indeed the transformation should fix this issue and get back the right range in this case!

I hope this helps!

[yellowcap]

Interesting! So what I thougt were bugs are really features. Thanks for the clarification @mathleur this makes sense now.

If you don't mind me asking, why are these conventions so different from the most common lat/lon grids? Latitude is almos always the Y dimension in mapping, and I would say most folks that are not from the weather world would expect longitude values to be in [-180, 180]. In other words, what is the coordinate reference system that this data is using?

[mathleur]

For the lat/lon grids, I believe we are currently following the ISO 6709 convention, where latitude is specified before longitude, but we are working towards changing the axes to latitude and longitude to make this clearer.

For the longitude ranges, we are following the GRIB format.

The data for the Destination Earth Climate Digital Twin is on a HEALPix grid and more information about the parameters and model can be found here.

I hope this helps answer the questions!

[yellowcap]

Yes thanks for all your help @mathleur and @awarde96 , it is very much appreciated.

The question of lat first vs lon first has been a source of counfusion for a while and is more to be expected.

However, the x and y mapping does not follow th ISO convention you mentioned. The wikipedia article reads

First horizontal position coordinate (Φ or y), such as latitude

Second horizontal coordinate (λ or x), such as longitude

So also here, y and x are lat / lon.

I'll close this ticket given that you are refactoring the coordinate handling anyway. But I would strongly suggest to map latitude to Y and longitude to X in the future.