SciTools/iris

Mapping a vector quantity from a curvilinear to rectilinear grid

DamienIrving opened this issue · 5 comments

DamienIrving commented

I initially posted the following as a question on the iris Google Group but figure it should actually be a GitHub issue since it triggered no error or warning when I ran it...

I'm currently working with CMIP5 ocean heat transport data. In particular, I need to calculate the northward heat transport at each latitude using the hfy (ocean_heat_y_transport) and hfx (ocean_heat_x_transport) variables. For most of the models, the ocean data is on a curvilinear grid, e.g:

print(hfy_cube)

ocean_heat_y_transport                     (time: 1872; cell index along second dimension: 384; cell index along first dimension: 320)
     Dimension coordinates:
          time                                  x                                        -                                      -
          cell index along second dimension     -                                        x                                      -
          cell index along first dimension      -                                        -                                      x
     Auxiliary coordinates:
          latitude                              -                                        x                                      x
          longitude                             -                                        x                                      x

In order to get the hfy data on a regular lat/lon grid, I was thinking I could use the following iris functionality,

target_grid = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
new_hfx_cube, new_hfy_cube = iris.analysis.cartography.rotate_winds(hfx_cube, hfy_cube, target_grid)

and then regrid to a rectilinear grid using iris.experimental.regrid.regrid_weighted_curvilinear_to_rectilinear.

When I go ahead and use iris.analysis.cartography.rotate_winds in this manner the values of new_hfx_cube remain unchanged from hfx_cube. I'm assuming this is because the original hfx and hfy cubes (created directly from a netCDF CMIP5 data file) have no coord_system:

print(hfx_cube.coord_system())
None

Does anyone know (a) if this lack of a coordinate system is the problem, and (b) if so, how do I specify a curvilinear grid? It doesn't appear to be an iris.coord_systems option...
http://scitools.org.uk/iris/docs/v1.9.2/iris/iris/coord_systems.html

corinnebosley commented

Hi Damien,
Have you tried checking the coordinate system on a specific coordinate of the cube rather than the whole thing?

Having had a little look at the iris functions you are using, it sounds like this might be something worth sending to the AVD Support Desk (ml-avd-support@metoffice.gov.uk) with a sample data file so that we can dig into it properly. Would you be able to do that?

pelson commented

Yes, rotate_winds assumes the x & y are lon/lats unless there is a coordinate system:

https://github.com/SciTools/iris/blob/master/lib/iris/analysis/cartography.py#L1035

if x_coord.coord_system is not None:
        src_crs = x_coord.coord_system.as_cartopy_projection()
    else:
        src_crs = ccrs.PlateCarree()

Out of interest, what is the coordinate system? Is it an ORCA grid?

FWIW cell index along second dimension is usually called projection_x_coordinate / projection_y_coordinate in CF terms (at least, that is how Iris has interpreted that part of the conventions).

pp-mo commented

Is it an ORCA grid?

We have been doing some work on regridding ORCA winds data internally here, and some progress has been made, but there are no plans to be include it in Iris at this point.

For Iris, there are specific problems with correctly handling this data, because of the specifics of the ORCA grids :

  • Although they aren't entirely "unstructured", they are defined only by 2d lat and lon arrays
  • The 2d structure implies cell connectivity, but it is not complete in that there are discontinuities (which in practice are always covered by missing data points)
  • The 2d arrays of lats and lons are definitive : there is no equational form to translate between geolocations and grid locations ...
  • ... accordingly, no continuous, separable 'native' X and Y coordinates' can be defined (unlike data on a map projection).

From Iris point of view there are also problems with the generality or specificity of operations to handle this : The connectivity of these grids is not explicitly recorded in the datafiles, and is effectively known only from provenance. However, this is necessary information for analysis operations, so we will need to invent a way for the user to specify that the data is of this type. Likewise, the exact grid on which winds are defined may also need to be asserted : I am told that NEMO and CICE data use different Arakawa grids with different cell locations of U and V values, and again this information is not explicit in the datafiles.

I hope this casts some light on why we have been slow to support more operations on this type of data !

github-actions commented

In order to maintain a backlog of relevant issues, we automatically label them as stale after 500 days of inactivity.

If this issue is still important to you, then please comment on this issue and the stale label will be removed.

Otherwise this issue will be automatically closed in 28 days time.

github-actions commented

This stale issue has been automatically closed due to a lack of community activity.

If you still care about this issue, then please either:

  • Re-open this issue, if you have sufficient permissions, or
  • Add a comment pinging @SciTools/iris-devs who will re-open on your behalf.