Create geometry for global regular lonlat from -180 to 180

[weiwilliam]

Hello,

I have a question about how to create the geometry with -180 to 180 in longitude.
I currently setup the geometry by following keys,

  function space: StructuredColumns
  grid:
    type: regular_lonlat
    nx: 576
    ny: 361
    domain :
      type : "rectangular"
      units : "degrees"
      ymax : 90.
      xmin : -180.
      ymin : -90.
      xmax : 180.
  halo: 1

It does give me the -180 to 180 geometry, however, the code below we used to prevent using missing value near the boundary seems not working.

atlas::mesh::actions::build_nodes_parallel_fields(mesh_);
atlas::mesh::actions::build_periodic_boundaries(mesh_);
atlas::mesh::actions::build_halo(mesh_, 1);

I wonder if there is other way to create the geometry or any solution to the missing value.

Thank you!

[wdeconinck]

Dear @weiwilliam without more information I cannot fully understand the issue.
What do you mean with "geometry"? How is the "mesh" created? What do you mean with "missing values"?
Could you provide a working code example with what you want to achieve?

[weiwilliam]

@wdeconinck Sorry about the incomplete information.
I collected more information yesterday, here are more details.
The geometry I referred to is the grid of regular lonlat based on the setting.

I want to access the model data stored in [-180, 180] and get the information at the observation location.
Originally, I only set up the grid with

    type: regular_lonlat
    nx: 576
    ny: 361

It creates the grid [0, 360], so the model data at the location is shifted 180 degrees.

So, I found adding the setting below can provide the grid properly.

    domain :
      type : "rectangular"
      units : "degrees"
      ymax : 90.
      xmin : -180.
      ymin : -90.
      xmax : 180.

But it has an issue when the observation is beyond the last grid point (179.325 in my case).

I output the grid lonlat information from mesh and functionspace like the way in issue #200, but we do add lines below in our code.

atlas::mesh::actions::build_nodes_parallel_fields(mesh_);
atlas::mesh::actions::build_periodic_boundaries(mesh_);
atlas::mesh::actions::build_halo(mesh_, 1);

For the grid [0, 360], it shows ghost points in the mesh_points.dat.
But it is not the case for grid [-180, 180], which only has 207,936 grid points for my dataset in mesh_points.dat.

I'm preparing the example code, and I will upload it later.

An alternative way,
atlas-meshgen 0to360.yaml -o mesh_0to360 --halo=1 --ghost --info
atlas-meshgen negpos180.yaml -o mesh_180 --halo=1 --ghost --info
Then you can see the differences between the ghost points section in mesh_0to360_info.msh and mesh_180_info.msh.
The prior has 209019 points, but the latter only has 207936 points and no grid points were marked as ghost.

Please let me know if you need more information.
Thank you.

[fmahebert]

Hi @wdeconinck! Just to clarify the above: @weiwilliam is working with JEDI and trying to adapt the atlas grid (in this case we're building an atlas lonlat grid and setting up a JEDI geometry wrapper on top of it) to match a dataset he has access to that is stored with longitudes ranging from [-180,180].

What Shih-Wei is finding is that when shifting the atlas lonlat grid from its default range of [0,360] to the new range of [-180,180] (by using the domain in Shih-Wei's yaml snippet above), there seems to be inconsistency in how the halos are set up. When the functionspace is NodeColumns, it appears the halos are added as expected and the JEDI interpolations are working; when the functionspace is StructuredColumns, it appears the halos are missing and JEDI's interpolations fail.

I realize the above is not yet very specific or actionable. If it's helpful, I can try to dig a bit to provide more specific details and a reproducer.

[wdeconinck]

Thanks for the clarification. An alternative snippet for the domain, if the domain is global after all, would be:

            domain :
                type : "global"
                units : "degrees"
                west : -180.

[weiwilliam]

Thank you, @wdeconinck!
With the snippet, JEDI interpolation is working with StructuredColumns.

[wdeconinck]

Interesting; even though I would recommend to use the "global" domain type, I'd still consider this a bug that needs fixing.
I created an issue to track this: #282

[wdeconinck]

And PR #283 should fix the problem, so that the previous code should work; though still just to change the 'west', I advise to use the global domain.

[fmahebert]

Thanks @wdeconinck for the advice and the quick bugfix!