Periodic grids

[matthew...@metoffice.gov.uk]

As part of the LFRic project <http://www.ecmwf.int/sites/default/files/elibrary/2016/16815-introduction-lfric-project.pdf> we are using UGRID to describe the surface grid of our model. This grid is then extruded to obtain a 3D mesh.

In order to avoid poles we have chosen the cubesphere for whole-planet models, however our scientists often find it useful to be able to test on a biperiodic Cartesian grid. Describing the periodic topology is not a problem with UGRID but we are not sure how to handle the co-ordinate space around the wrap-around.

The width of a cell is determined from the co-ordinates at each end but this breaks down at the edges. Here that calculation goes massively negative. This leaves us with the situation of not knowing how large the final cell is. We need this information as the grid may not be regular.

We could store the width of each cell as a separate data entity but this leads to large duplication. We would be storing information we can infer for all but the border cells. Ideally we would store the width information only for those border cells but I'm not sure if this is supported by UGRID at present.

I have to assume that people are using periodic grids all over the place and that this problem is already solved. So if someone could point me in the right direction I would be grateful.

--
Matthew Hambley
  Met Office Hadley Centre, Exeter, UK

[David Ham]

Dear Matthew,

There are a few approaches which are frequently employed for the representation of periodic domains on unstructured meshes.

1. Represent the coordinate field in a discontinuous function space. This way each vertex can have different coordinates when seen from different cells. Note that LFRic is already committed to supporting discontinuous spaces for tracers, so discontinuous coordinates is not a big stretch. The key advantage of this approach is that it removes all special casing of periodic domains. It does require the storage of more data than other approaches, however this additional storage is likely to be small in comparison with time-varying quantities. With reference to Gung Ho and LFRic, you might be interested to know that this is the approach which Firedrake takes.

2.  Store two topologies, one periodic and one not. Store the coordinates on the non-periodic topology and the solution fields on the periodic topology. This also avoids explicit special casing in the dynamics code itself, but does require some special purpose bookkeeping code in your LFRic setup.  It does store less data than option 1. This is the approach that Fluidity takes (http://fluidityproject.github.io).

3. Explicitly store a flagged boundary on your cells and put special purpose code in your dynamics to change coordinates when you hit the problematic cells. This is rather intrusive as it requires special purpose code in the dynamics, although if you really want to go this route, it is plausible that this code could be included in LFRic's code generation layer.

Coincidentally, I happen to be at the Met Office today. If you would like to meet in person for a higher bandwidth discussion of the various ways to do this, I'm very happy to do so.

Regards,

David

--
You received this message because you are subscribed to the Google Groups "UGRID Interoperability" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ugrid-interoperab...@googlegroups.com.
To post to this group, send email to ugrid-inter...@googlegroups.com.
Visit this group at https://groups.google.com/group/ugrid-interoperability.
For more options, visit https://groups.google.com/d/optout.

--

Dr David Ham

Department of Mathematics

Imperial College London

[Hambley, Matthew]

From: David Ham [mailto:da...@ham.dropbear.id.au]
[snip]

> There are a few approaches which are frequently employed for the
> representation of periodic domains on unstructured meshes.

Yes, that is how these things are represented in memory. But I am asking how to best represent them in a UGRID file.

--
Matthew Hambley
Scientific Software Engineer

[Chris Barker]

It seems to me that most of the periodic issues are with how you interpret and process the grid, not how it is represented in a file. So what am I missing? what information do you want to put in the file that you can't with the existing convention?

i.e. the face_face_connectivity array specified which faces are neighbors of which faces, to a face can be the neighbor of a face that's on the "other side" of the world just fine.

 we are not sure how to handle the co-ordinate space around the wrap-around.

The width of a cell is determined from the co-ordinates at each end but this breaks down at the edges. Here that calculation goes massively negative.

IIUC, this is the classic GIS problem with wrap around of polygons (or polylines) that cross the dateline. -- a line from -179 to 179 -- is that 2 degrees or 358 degrees long? But that seems like a problem that needs to be solved in your model/processing/visualizing software?

i.e. your code has a simple check -- if the width calculation goes negative, then you do the math "the other way around"

But yes, storing that info in the file would be nice.

We could store the width of each cell as a separate data entity but this leads to large duplication. 

We would be storing information we can infer for all but the border cells. Ideally we would store the width information only for those border cells but I'm not sure if this is supported by UGRID at present.

I didn't think so either. But I just took a look, and found "location index set":

http://ugrid-conventions.github.io/ugrid-conventions/#location-index-set

Could you use this to specify the cell width for only those cells where it is not obvious?

-CHB

--

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris....@noaa.gov

[Hambley, Matthew]

From: ugrid-inter...@googlegroups.com [mailto:ugrid-

>
> It seems to me that most of the periodic issues are with how you
> interpret and process the grid, not how it is represented in a file. So
> what am I missing? what information do you want to put in the file that
> you can't with the existing convention?

Let me see if I can explain more clearly.

For each cell in the grid we can calculate the edge length (lets assume a 1D grid for simplicity) from the co-ordinates of the node at each end. Except, this can't be done for the edge between the last node and the first node, i.e. the wrap-around.

Thus we are left not knowing how long the wrap-around edge is.

The question is, how do we utilise a UGRID file so this information is available. There is no obvious solution so I was hoping the list could point me in the direction of a nonobvious solution.

--
Matthew Hambley
Scientific Software Engineer

[Christopher Barker]

On Wed, Dec 14, 2016 at 12:38 AM Hambley, Matthew <matthew...@metoffice.gov.uk> wrote:

>For each cell in the grid we can calculate the edge length (lets assume a 1D grid for simplicity) from the co-ordinates of the node at each end. Except, this can't be done for the edge between the last node and the first node, i.e. the wrap-around.



> Thus we are left not knowing how long the wrap-around edge is.

Got it. I've delt with this in program logic. On a globe, there are two great circle distances between two points. If you can assume the one you want is the shorter one, which it would invariably be in this case, then you have ugly logic, but the result is well defined.

>The question is, how do we utilise a UGRID file so this information is available. 

So you may want to avoid the expensive logic in your program code. And when you write the file, you have extra information.

So I'll put it back on you: how would you express this information in ANY file format or convention, or in your processing code.

For instance, I think this is an unsolved problem in the various GIS conventions.

But I can think of a couple options:

1) You store a flag on the edge, or face, or ... indicating that it crosses the wrap-around. 

2) you store the size of the edges, faces etc. so they don't need to be computed.

I'd probably do that as ordinary "data on the grid" variables -- probably not that big compared to other data, and a flag at least would compress really well.

But if you want to only store info on the parts of the grid that overlap the wrap around, you can do that with a "location index set".

As there are multiple ways to handle this, perhaps we should establish a recommended way in the convention.

I don't generally work with global scale models, to I don't have an opinion.

But I'd say we start with looking at other conventions-- this isn't a UGRID specific issue. Is this addressed in CF at all? In any of the osgeo conventions?

HTH,

-CHB

[Hambley, Matthew]

From: Christopher Barker [mailto:pyth...@gmail.com]

>
>>For each cell in the grid we can calculate the edge length
>> (lets assume a 1D grid for simplicity) from the co-ordinates of the
>> node at each end. Except, this can't be done for the edge between the
>> last node and the first node, i.e. the wrap-around.
>>
>> Thus we are left not knowing how long the wrap-around edge is.
>
> Got it. I've delt with this in program logic. On a globe,
> there are two great circle distances between two points. If you can
> assume the one you want is the shorter one, which it would invariably
> be in this case, then you have ugly logic, but the result is well
> defined.

That is true for a globe where you can assume 360degrees around a circle. However we are dealing with a biperiodic Cartesian domain. In other words it is size 'n' by 'm' which wraps around at all edges. It is the surface of a 3D torus but has 2D coordinates, if that makes sense.

>>The question is, how do we utilise a UGRID file so this
>> information is available.
>
> So you may want to avoid the expensive logic in your program
> code. And when you write the file, you have extra information.

That is the case. Our mesh generator knows how big the wrap around cells are. We just need to communicate that information to our model via the UGRID file.

[snip]

> For instance, I think this is an unsolved problem in the
> various GIS conventions.

If this is true then we can invent something. I just find it hard to believe that we are the only group in the world using this kind of grid.

If I tackle your two suggestions in reverse order:

> 2) you store the size of the edges, faces etc. so they don't
> need to be computed.

This is the obvious solution but we don't like this because it means that for most nodes you have two sources of truth regarding the length of an edge. You have the stored value and the value computed from end point coordinates. We would rather avoid this duplication so as to avoid the possibility of them disagreeing and for reasons of file bloat.

> 1) You store a flag on the edge, or face, or ... indicating
> that it crosses the wrap-around.

This, then, is the next obvious option. The question is how to achieve it.

The UGRID standard document talks about "boundary_node_connectivity" and the ability to tag this boundary as either "open" or "closed". Do either of these refer to wrap-around? If so it continues to talk about associating "edge_coordinates" with the boundary and how these coordinates may be associated with "bounds". If this describes what I think it does then those bounds would describe the length of the edge.

[snip]

> But if you want to only store info on the parts of the grid
> that overlap the wrap around, you can do that with a "location index
> set".

This looks like it might be simpler than what I describe above. However the "boundary_node_connectivity" approach looks like it might be the method which was intended for this purpose. We want to avoid misappropriating things in non-standard ways.

> As there are multiple ways to handle this, perhaps we should
> establish a recommended way in the convention.

Given how common biperiodic Cartesian domains are it seems like a worthwhile thing to standardise. After all "Spacewar!" was doing this on a PDP-1 in 1962. :-)

[snip]

> But I'd say we start with looking at other conventions-- this
> isn't a UGRID specific issue. Is this addressed in CF at all? In any of
> the osgeo conventions?

I assume you are referring to the NetCDF climate and forecast metadata convention. In which case, no, they seem to deal exclusively with latitude/longitude. I had a quick look at the Open Source Geospatial Foundation's web page but couldn't find any obvious conventions to look at.

[Chris Barker]

On Thu, Dec 15, 2016 at 3:04 AM, Hambley, Matthew <matthew...@metoffice.gov.uk> wrote:

> If you can
> assume the one you want is the shorter one, which it would invariably
> be in this case, then you have ugly logic, but the result is well
> defined.

That is true for a globe where you can assume 360degrees around a circle. However we are dealing with a biperiodic Cartesian domain. In other words it is size 'n' by 'm' which wraps around at all edges. It is the surface of a 3D torus but has 2D coordinates, if that makes sense.

sure -- but the principle is the same -- there are two "directions" you can compute a distance over -- it's not a 360 degrees, but a cell is presumable much smaller than the domain, so it's deterministic how to compute the distance over the wrap-around. nonethe less, you may not want to have to do that.

>>The question is, how do we utilise a UGRID file so this
>> information is available.
>
>               So you may want to avoid the expensive logic in your program
> code. And when you write the file, you have extra information.

That is the case. Our mesh generator knows how big the wrap around cells are. We just need to communicate that information to our model via the UGRID file.

[snip]
>               For instance, I think this is an unsolved problem in the
> various GIS conventions.

If this is true then we can invent something. I just find it hard to believe that we are the only group in the world using this kind of grid.

I'm sure not -- and other global grids deal with this one way or another -- anyone?

If I tackle your two suggestions in reverse order:

>               2) you store the size of the edges, faces etc. so they don't
> need to be computed.

This is the obvious solution but we don't like this because it means that for most nodes you have two sources of truth regarding the length of an edge. You have the stored value and the value computed from end point coordinates. We would rather avoid this duplication so as to avoid the possibility of them disagreeing and for reasons of file bloat.

yeah -- I sympathize, though UGRID does often store redundant information -- a lot of the connectivity arrays can be computed on the fly -- as well as the coordinate arrays. So it's not without precedent.

and if you store it all, then it's easy to extract the info -- rather than having to have a check in there:

if this_is_a_wrapped_cell:

    read_size_from_file

else:

    compute_size_from_coordinates 

>               1) You store a flag on the edge, or face, or ... indicating
> that it crosses the wrap-around.

This, then, is the next obvious option. The question is how to achieve it.

The UGRID standard document talks about "boundary_node_connectivity" and the ability to tag this boundary as either "open" or "closed". Do either of these refer to wrap-around? If so it continues to talk about associating "edge_coordinates" with the boundary and how these coordinates may be associated with "bounds". If this describes what I think it does then those bounds would describe the length of the edge.

the idea behind boundary_node_connectivity is to be able to define a subset of edges -- typically the boundary of the mesh -- so that you can then specify data on those edges -- I don't think the standard restricts what data is set on those boundaries, or how to specify them.

from the doc:

"""
Information about the nature of each boundary edge (e.g. open/closed, land/water, grouping, etc.) may optionally be stored in ancillary boundary-type variables
"""

note the open/closed is an "e.g." -- just an example -- you can put any data you want on a  variable (you should follow CF conventions for that variable, though.

so if you want to specify a flag for a wrapped boundary, this would be a good way to do it.

And you could have an edge_length variable with the length, tied only to the edges specified by the boundary_node_connectivity (maybe we shouldn't have used the term boundary -- "edge_subset"???

[snip]
>               But if you want to only store info on the parts of the grid
> that overlap the wrap around, you can do that with a "location index
> set".

This looks like it might be simpler than what I describe above. However the "boundary_node_connectivity" approach looks like it might be the method which was intended for this purpose. We want to avoid misappropriating things in non-standard ways.

location_index_set  allows you to specify a subset of some other complete set of locations -- subset of nodes, edges, faces, etc.

boundaries are only edges (hmm in 3d, you might want faces as boundaries...) and may be separately specified without having to have the full set of edges defined -- in common cases, the nodes and faces are defined, and there is no need for a full set of edges -- redundant with the cell info -- yes?

So if you have a complete edge_node_connectivity array, then you could use location_index_set, but if not, you would need to use boundary_node_connectivity.

>               As there are multiple ways to handle this, perhaps we should
> establish a recommended way in the convention.

Given how common biperiodic Cartesian domains are it seems like a worthwhile thing to standardise.

yup -- once you hae thi figured out -- maybe propose an addition to the standard? Ideally, we'd get a few other folks with similar models to comment.

 

I assume you are referring to the NetCDF climate and forecast metadata convention. In which case, no, they seem to deal exclusively with latitude/longitude.

no it's not always lat=-ln -- other projections can be used -- and even with lat-lon, you have similar wrap-around issues on the globe.

but yes -- I haven't seen anything about dealing with this.

I had a quick look at the Open Source Geospatial Foundation's web page but couldn't find any obvious conventions to look at.

yeah -- I haven't seen anything -- it seems the GIS community mostly "punts" by shifting coordinates (or choosing a projection) where the "wrap around" is outside the domain of interest -- but you can't do that with a global dataset.

In our case, we set up our oil spill model to handle longitude from -180 to + 360 degrees. Then it's up to the user to get all the input data in the right system -- we have no idea that -180 and +180 are the same point.

and it seems that Web Mapping, and everything else I've seen does something like that.

but that simply wouldn't work for something moving all the way around the earth...

Hmm -- I wonder what the airplane route folks do?

-CHB

[Hambley, Matthew]

From: Chris Barker [mailto:chris....@noaa.gov]

> On Thu, Dec 15, 2016 at 3:04 AM, Hambley, Matthew
> <mailto:matthew...@metoffice.gov.uk> > wrote:
>>
>>> If you can assume the one you want is the shorter one, which it would
>>> invariably be in this case, then you have ugly logic, but the result is
>>> well defined.
>>
>> That is true for a globe where you can assume 360degrees around a
>> circle. However we are dealing with a biperiodic Cartesian domain. In
>> other words it is size 'n' by 'm' which wraps around at all edges. It
>> is the surface of a 3D torus but has 2D coordinates, if that makes
>> sense.
>
> sure -- but the principle is the same -- there are two "directions" you
> can compute a distance over -- it's not a 360 degrees, but a cell is
> presumable much smaller than the domain, so it's deterministic how to
> compute the distance over the wrap-around. nonethe less, you may not
> want to have to do that.

That would be true if we knew the size of the domain and we were always dealing with rectangular meshes. However the point is that we don't know the size of the domain, we wish to determine it from the UGRID file. Plus we may well want to play with 2D triangular or hexagonal meshes in the future.

[snip]

> and if you store it all, then it's easy to extract the info -- rather
> than having to have a check in there:
>
> if this_is_a_wrapped_cell:
> read_size_from_file
> else:
> compute_size_from_coordinates

While we look to avoid conditional statements in our time step code for performance reasons, we don't really care about them in initialisation.

[snip]

>>> As there are multiple ways to handle this, perhaps we should
>>> establish a recommended way in the convention.
>>
>> Given how common biperiodic Cartesian domains are it seems like a
>> worthwhile thing to standardise.
>
> yup -- once you hae thi figured out -- maybe propose an addition to the
> standard? Ideally, we'd get a few other folks with similar models to
> comment.

After further discussion with colleagues we have come to the conclusion that the most flexible approach, which works for all possible meshes (that we can think of), is to have a location index set of duplicate vertices. These will have two pieces of information: Co-ordinates and an index into the array of vertices saying "that's me over there".

In this scheme there are duplicates of the "left hand" vertices at the "right hand" of the mesh. Topologically these are the same vertex but geometrically they represent the "n+1" version of it.

Does this fit with how UGRID is envisaged to work? Does it make sense? Is there a better way?

--
Matthew Hambley - Dark Lord of the LFRic build system

  Scientific Software Engineer
  Met Office Hadley Centre, Exeter, UK

...I have a ticket for that

[Chris Barker]

On Fri, Jan 6, 2017 at 3:27 AM, Hambley, Matthew <matthew...@metoffice.gov.uk> wrote:

> sure -- but the principle is the same -- there are two "directions" you
> can compute a distance over -- it's not a 360 degrees, but a cell is
> presumable much smaller than the domain, so it's deterministic how to
> compute the distance over the wrap-around. nonethe less, you may not
> want to have to do that.

That would be true if we knew the size of the domain and we were always dealing with rectangular meshes. However the point is that we don't know the size of the domain, we wish to determine it from the UGRID file. Plus we may well want to play with 2D triangular or hexagonal meshes in the future.

I'm still confused -- but you know your use case.

 

[snip]
> and if you store it all, then it's easy to extract the info -- rather
> than having to have a check in there:
>
> if this_is_a_wrapped_cell:
>     read_size_from_file
> else:
>     compute_size_from_coordinates

While we look to avoid conditional statements in our time step code for performance reasons, we don't really care about them in initialisation.

I think it's less a performance issue than a code complexity one - if these details have already been worked out, it's nice to have the information in the file.

After further discussion with colleagues we have come to the conclusion that the most flexible approach, which works for all possible meshes (that we can think of), is to have a location index set of duplicate vertices. These will have two pieces of information: Co-ordinates and an index into the array of vertices saying "that's me over there".

In this scheme there are duplicates of the "left hand" vertices at the "right hand" of the mesh. Topologically these are the same vertex but geometrically they represent the "n+1" version of it.

Does this fit with how UGRID is envisaged to work? Does it make sense? Is there a better way?

This sounds good to me.

I suggest you post a (ideally small) example file or CDL version, to make sure we're all clear on what it means.

And then we can add it to the Doc as an example.

-Chris