Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Can't get area weights of cube with additional auxillary coordinates of lat and lon #3916

Closed
nhsavage opened this issue Nov 3, 2020 · 3 comments · Fixed by #5029
Closed

Comments

@nhsavage
Copy link
Contributor

nhsavage commented Nov 3, 2020

🐛 Bug Report

Iris is unable to calculate area weights for rotated pole grids if auxiliary coordinates with true latitude and longitude have been added to the cube.

How To Reproduce

  1. load data on a rotated pole grid
# load rotated pole grids
fname = iris.sample_data_path("rotated_pole.nc")
cube = iris.load_cube(fname)
  1. attach the aux coords thus:
# get cube's coordinate system
cs = cube.coord_system()
# get coord names
# Longitude
xcoord = cube.coord(axis="X", dim_coords=True)
# Latitude
ycoord = cube.coord(axis="Y", dim_coords=True)

# read in the grid lat/lon points from the cube
glat = cube.coord(ycoord).points
glon = cube.coord(xcoord).points

# create a rectangular grid out of an array of
# glon and glat values, shape will be len(glat)xlen(glon)
x, y = np.meshgrid(glon, glat)

# get the cube dimensions which corresponds to glon and glat
x_dim = cube.coord_dims(xcoord)[0]
y_dim = cube.coord_dims(ycoord)[0]

# define two new variables to hold the unrotated coordinates
rlongitude, rlatitude = iris.analysis.cartography.unrotate_pole(
    x, y, cs.grid_north_pole_longitude, cs.grid_north_pole_latitude
)

# create two new auxillary coordinates to hold
# the values of the unrotated coordinates
reg_long = iris.coords.AuxCoord(rlongitude, long_name="longitude", units="degrees")
reg_lat = iris.coords.AuxCoord(rlatitude, long_name="latitude", units="degrees")

# add two auxilary coordinates to the cube holding
# regular(unrotated) lat/lon values
cube.add_aux_coord(reg_long, [y_dim, x_dim])
cube.add_aux_coord(reg_lat, [y_dim, x_dim])

  1. print the cube:
air_temperature / (K)               (grid_latitude: 550; grid_longitude: 700)
     Dimension coordinates:
          grid_latitude                           x                    -
          grid_longitude                          -                    x
     Auxiliary coordinates:
          latitude                                x                    x
          longitude                               x                    x

4. try and get area weights:
weights =iris.analysis.cartography.area_weights(cube)

Traceback (most recent call last):
File "area_bug.py", line 49, in
weights =iris.analysis.cartography.area_weights(cube)
File ".../lib/python3.6/site-packages/iris/analysis/cartography.py", line 399, in area_weights
lon, lat = _get_lon_lat_coords(cube)
File ".../python3.6/site-packages/iris/analysis/cartography.py", line 189, in _get_lon_lat_coords
"Calling _get_lon_lat_coords with multiple lat or lon coords"
ValueError: Calling _get_lon_lat_coords with multiple lat or lon coords is currently disallowed

5. remove the extra coords - works fine:

cube.remove_coord('longitude')
cube.remove_coord('latitude')
weights =iris.analysis.cartography.area_weights(cube)

Expected behaviour

I expect that a cube like this should return the correct area weights

Environment

  • RHEL 7.9
  • Iris Version: 2.2.0

Additional context

Click to expand this section...
https://github.com/SciTools/iris/blob/master/lib/iris/analys... 

I would suggest that this: 

coord for coord in cube.coords() if "latitude" in coord.name() 

 

should actually be: 

coord for coord in cube.coords(dim_coords=True) if "latitude" in coord.name()

See here for further details.

@pp-mo
Copy link
Member

pp-mo commented Nov 11, 2020

Don't want to overcomplicate, but this reminds me of : #3169
In an ideal world... instead of returning an array of values, the calculation could add a cell-measure to the cube, and the stats methods could use that automatically. So, I wish we could do that.

A separate issue really.
Just saying..

@trexfeathers trexfeathers added this to To do in Peloton via automation Jan 6, 2021
@trexfeathers trexfeathers added the Peloton 🚴‍♂️ Target a breakaway issue to be caught and closed by the peloton label Jan 6, 2021
@alastair-gemmell
Copy link
Contributor

Hi @nhsavage, just checking through some older issues and wonder to what degree this is still relevant to you?

@nhsavage
Copy link
Contributor Author

Hi @nhsavage, just checking through some older issues and wonder to what degree this is still relevant to you?
@alastair-gemmell - we do have work arounds for this and it's generally not an issue. It is a common with roated pole data to add these auxilliary coordinates as it makes it easier to work with for many users. On balance - not that important but a "nice to have"

@trexfeathers trexfeathers removed the Peloton 🚴‍♂️ Target a breakaway issue to be caught and closed by the peloton label Apr 6, 2022
@trexfeathers trexfeathers removed this from Backlog in Peloton Apr 6, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
Development

Successfully merging a pull request may close this issue.

5 participants