Re-grid cube to a defined grid

[Patrick Breach]

I know I can regrid one cube to that of another via,

cube1.regrid(cube2, iris.analysis.Linear())

But how can I define a grid that is not necessarily connected with another cube? Say I have a  cube ("cube1"), with a 4 x 4 degree rectilinear grid and want to regrid to 2.5 x 2.5 degrees

I am going with this approach here so far,

import numpy as np
import pandas as pd
from iris.pandas import as_cube

#Use the extends of one cube to generate a set of lat/lon coordinate dimensios
lat,lon = cube1.coord('latitude'),cube1.coord('longitude')
lat_min,lon_min = lat.points.min(),lon.points.min()
lat_max,lon_max = lat.points.max(),lon.points.max()

lat25 = np.arange(lat_min, lat_max, 2.5)
lon25 = np.arange(lon_min, lon_max, 2.5)

#Generate a 'grid cube' from a pandas DataFrame to act as the target grid
df = pd.DataFrame(np.random.random((len(lat25),len(lon25))))
df.index = lat25
df.columns = lon25

grid_cube = as_cube(df)
grid_cube.coord('index').rename('latitude')
grid_cube.coord('columns').rename('longitude')
grid_cube.coord('latitude').units = 'degrees'
grid_cube.coord('longitude').units = 'degrees'

#Do the regridding
cube1.regrid(grid_cube, iris.analysis.Linear())

Which works, however this seems to be inefficient and extremely verbose. Is there an easier way to do this? If not shouldn't there be?

[Andrew Dawson]

Hi Patrick

I think in Iris the term regridding specifically refers to interpolating a cube onto the grid of another cube (or at least in the context of the regrid method it does). I think what you want is probably the interpolate method. You should be able to do something like this:

lat,lon = cube1.coord('latitude'),cube1.coord('longitude')
lat_min,lon_min = lat.points.min(),lon.points.min()
lat_max,lon_max = lat.points.max(),lon.points.max()
lat25 = np.arange(lat_min, lat_max, 2.5)
lon25 = np.arange(lon_min, lon_max, 2.5)

result = cube1.interpolate([('latitude', lat25), ('longitude', lon25)],
                           iris.analysis.Linear())

I hope that helps, if not let us know.

The User Guide chapter on interpolation and regridding may be a useful read: http://scitools.org.uk/iris/docs/latest/userguide/interpolation_and_regridding.html

[Patrick Breach]

Andrew,

This is much simpler! I read through the chapter previously but I guess I didn't think of using interpolate this way.

Thanks,

Patrick

[LeonH]

I'm trying this approach as regrettably setting the central_longitude keyword in my projection that I call iris.analysis.cartography.project() with has no effect, I still end up with a cube that has as its first longitude -179.5, when what I want is zero. No worry, maybe I can use interpolate! However, I must be doing something wrong as this function:

def startAtZero(cube):
    """
    The project function above unfortunately ignores its keyword
    so I have to interpolate the data manually
    so that it fits with what I want
    """
    lonPoints = cube.coord('longitude').points - 179.5
    cube.coord('longitude').circular = True
    samplePoints = [('longitude', lonPoints)]
    return cube.interpolate(samplePoints, iris.analysis.Linear())

Has this error:

>>> iris.__version__
'1.7.3-rc.1'
>>> t = op.startAtZero(test)
Traceback (most recent call last):
  File "<pyshell#94>", line 1, in <module>
    t = op.startAtZero(test)
  File "orca_project.py", line 53, in startAtZero
    return cube.interpolate(samplePoints, iris.analysis.Linear())
  File "/project/avd/iris/live/pre-release/iris/cube.py", line 3231, in interpolate
    interp = scheme.interpolator(self, coords)
  File "/project/avd/iris/live/pre-release/iris/analysis/__init__.py", line 1561, in interpolator
    extrapolation_mode=self.extrapolation_mode)
  File "/project/avd/iris/live/pre-release/iris/analysis/_interpolator.py", line 162, in __init__
    self._setup()
  File "/project/avd/iris/live/pre-release/iris/analysis/_interpolator.py", line 339, in _setup
    offset = ((coord_points.max() + coord_points.min() - modulus)
TypeError: unsupported operand type(s) for -: 'float' and 'NoneType'
>>> print test
integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content / (J/m2) (latitude: 180; longitude: 360)
     Dimension coordinates:
          latitude                                                                                x               -
          longitude                                                                               -               x
     Scalar coordinates:
          time: 2010-01-16 12:00:00, bound=(2010-01-01 00:00:00, 2010-02-01 00:00:00)
     Attributes:
          Conventions: CF-1.1
          NCO: 4.2.0
          interval_operation: 1350.0
          interval_write: -1.0
          online_operation: ave(only(x))
          production: An IPSL model
>>> print test.coord('longitude')
DimCoord(array([-179.5, -178.5, -177.5, -176.5, -175.5, -174.5, -173.5, -172.5,
       -171.5, -170.5, -169.5, -168.5, -167.5, -166.5, -165.5, -164.5,
         ...
        164.5,  165.5,  166.5,  167.5,  168.5,  169.5,  170.5,  171.5,
        172.5,  173.5,  174.5,  175.5,  176.5,  177.5,  178.5,  179.5]), standard_name='longitude', units=Unit('1'), circular=True)

Can anyone see what I am doing wrong?

[Andrew Dawson]

The traceback indicates that the modulus attribute of the coordinates units attribute (coord.units.modulus) is set to None. The code will assume a value of modulus=0 if there is no modulus attribute present, but if there is one (even if it is None) it will be honoured. I notice the units of your longitude coordinate are '1' (dimensionless) which seems wrong, they should really be degrees_east from looking at the values, or at least degrees.

It seems like the immediate solution would be to set the units attribute correctly before interpolation.

Failing that, perhaps you could tell us a bit more about the problem. Things like where the data are loaded from and how, what processing has been applied etc. would be useful.