CF CoordSystems (replacement pull)

docs/iris/example_code/graphics/custom_file_loading.py

@@ -164,10 +164,7 @@ def NAME_to_cube(filenames, callback):
             cube.add_aux_coord(time_coord)
 
             # build a coordinate system which can be referenced by latitude and longitude coordinates
-            lat_lon_coord_system = icoord_systems.LatLonCS( icoord_systems.SpheroidDatum("spherical", 6371229.0, flattening=0.0, units='m'), 
-                                                     icoord_systems.PrimeMeridian(label="Greenwich", value=0.0), 
-                                                     n_pole=icoord_systems.GeoPosition(90, 0), reference_longitude=0.0
-                                                    )
+            lat_lon_coord_system = icoord_systems.GeogCS(6371229)
 
             # build regular latitude and longitude coordinates which have bounds
             start = header['X grid origin'] + header['X grid resolution']

docs/iris/src/userguide/cube_statistics.rst

@@ -108,7 +108,7 @@ These areas can now be passed to the ``collapsed`` method as weights:
          Scalar coordinates:
               grid_latitude: Cell(point=1.5145501, bound=(0.14430022, 2.8848)) degrees
               grid_longitude: Cell(point=358.74948, bound=(357.49399, 360.00497)) degrees
-              surface_altitude: Cell(point=399.625, bound=(-14.0, 813.25)) m
+              surface_altitude: Cell(point=399.625, bound=(-14.000001, 813.25)) m
          Attributes:
               STASH: m01s00i004
               history: Mean of air_potential_temperature over grid_longitude, grid_latitude

lib/iris/analysis/calculus.py

@@ -412,7 +412,7 @@ def curl(i_cube, j_cube, k_cube=None, ignore=None, update_history=True):
 
     Return (i_cmpt_curl_cube, j_cmpt_curl_cube, k_cmpt_curl_cube)
 
-    The calculation of curl is dependent on the type of :func:`iris.coord_systems.HorizontalCS` in the cube:
+    The calculation of curl is dependent on the type of :func:`iris.coord_systems.CoordSystem` in the cube:
 
         Cartesian curl
 
@@ -440,16 +440,6 @@ def curl(i_cube, j_cube, k_cube=None, ignore=None, update_history=True):
         ignore = None
         warnings.warn('The ignore keyword to iris.analysis.calculus.curl is deprecated, ignoring is now done automatically.')
 
-    # get the radius of the earth
-    latlon_cs = i_cube.coord_system(iris.coord_systems.LatLonCS)
-    if latlon_cs and latlon_cs.datum.is_spherical():
-        r = latlon_cs.datum.semi_major_axis
-        r_unit = latlon_cs.datum.units
-    else:
-        r = iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
-        r_unit = iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS_UNIT
-
-
     # Get the vector quantity names (i.e. ['easterly', 'northerly', 'vertical'])
     vector_quantity_names, phenomenon_name = spatial_vectors_with_phenom_name(i_cube, j_cube, k_cube)
 
@@ -478,11 +468,11 @@ def curl(i_cube, j_cube, k_cube=None, ignore=None, update_history=True):
 
     y_dim = i_cube.coord_dims(y_coord)[0]
 
-    horiz_cs = i_cube.coord_system('HorizontalCS')
-    if horiz_cs is None:
-        raise ValueError('Could not get the horizontal CS of the cubes provided.')
+    horiz_cs = i_cube.coord_system('CoordSystem')
 
-    if horiz_cs.cs_type == iris.coord_systems.CARTESIAN_CS:
+    # Planar (non spherical) coords?
+    ellipsoidal = isinstance(horiz_cs, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS))
+    if not ellipsoidal:
 
         # TODO Implement some mechanism for conforming to a common grid
         dj_dx = _curl_differentiate(j_cube, x_coord)
@@ -525,17 +515,32 @@ def curl(i_cube, j_cube, k_cube=None, ignore=None, update_history=True):
 
         result = [i_cmpt, j_cmpt, k_cmpt]
 
-    elif horiz_cs.cs_type == iris.coord_systems.SPHERICAL_CS:
+    # Spherical coords (GeogCS or RotatedGeogCS).
+    else:
         # A_\phi = i ; A_\theta = j ; A_\r = k
         # theta = lat ; phi = long ;
         # r_cmpt = 1/ ( r * cos(lat) ) * ( d/dtheta ( i_cube * sin( lat ) ) - d_j_cube_dphi )
         # phi_cmpt = 1/r * ( d/dr (r * j_cube) - d_k_cube_dtheta)
         # theta_cmpt = 1/r * ( 1/cos(lat) * d_k_cube_dphi - d/dr (r * i_cube)
-        if not horiz_cs.datum.is_spherical():
-            raise NotImplementedError('Cannot take the curl over a non-spherical datum.')
-
         if y_coord.name() != 'latitude' or x_coord.name() != 'longitude':
             raise ValueError('Expecting latitude as the y coord and longitude as the x coord for spherical curl.')
+
+        # Get the radius of the earth - and check for sphericity
+        ellipsoid = horiz_cs
+        if isinstance(horiz_cs, iris.coord_systems.RotatedGeogCS):
+            ellipsoid = horiz_cs.ellipsoid
+        if ellipsoid:
+            # TODO: Add a test for this
+            r = ellipsoid.semi_major_axis
+            r_unit = iris.unit.Unit("m")
+            spherical = (ellipsoid.inverse_flattening == 0.0)
+        else:
+            r = iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
+            r_unit = iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS_UNIT
+            spherical = True
+
+        if not spherical:
+            raise ValueError("Cannot take the curl over a non-spherical ellipsoid.")
 
         lon_coord = x_coord.unit_converted('radians')
         lat_coord = y_coord.unit_converted('radians')
@@ -589,11 +594,7 @@ def curl(i_cube, j_cube, k_cube=None, ignore=None, update_history=True):
         d_k_cube_dphi = dri_dr = None
 
         result = [phi_cmpt, theta_cmpt, r_cmpt]
-
-    else:
-        raise ValueError("Horizontal coord system neither cartesian nor spherical spheroid: %s %s (%s)" \
-                         % (type(horiz_cs), horiz_cs.cs_type, horiz_cs.datum))
-
+
     for direction, cube in zip(vector_quantity_names, result):
         if cube is not None:
             cube.rename('%s curl of %s' % (direction, phenomenon_name))

lib/iris/analysis/cartography.py

@@ -32,9 +32,10 @@
 import iris.unit
 
 
-#This value is used as a fall-back if the cube does not define the earth 
-DEFAULT_SPHERICAL_EARTH_RADIUS = 6367.47
-DEFAULT_SPHERICAL_EARTH_RADIUS_UNIT = iris.unit.Unit('kilometres')
+# This value is used as a fall-back if the cube does not define the earth 
+DEFAULT_SPHERICAL_EARTH_RADIUS = 6367470
+# TODO: This should not be necessary, as CF is always in meters
+DEFAULT_SPHERICAL_EARTH_RADIUS_UNIT = iris.unit.Unit('m')
 
 
 def wrap_lons(lons, base, period):
@@ -108,9 +109,14 @@ def lat_lon_range(cube, mode=None):
     used in the min/max calculation. (Must be one of iris.coords.POINT_MODE or iris.coords.BOUND_MODE)
 
     """
+    # Helpful error if we have an inappropriate CoordSystem
+    cs = cube.coord_system("CoordSystem")
+    if cs is not None and not isinstance(cs, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)):
+        raise ValueError("Latlon coords cannot be found with {0}.".format(type(cs)))
+
     # get the lat and lon coords (might have "grid_" at the start of the name, if rotated).
     lat_coord, lon_coord = _get_lat_lon_coords(cube)
-    cs = cube.coord_system('LatLonCS')
+    cs = cube.coord_system('CoordSystem')
 
     if lon_coord.has_bounds() != lat_coord.has_bounds():
         raise ValueError('Cannot get the range of the latitude and longitude coordinates if they do '
@@ -123,7 +129,7 @@ def lat_lon_range(cube, mode=None):
     else:
         _mode = iris.coords.POINT_MODE
 
-    if cs.has_rotated_pole():
+    if isinstance(cs, iris.coord_systems.RotatedGeogCS):
         if _mode == iris.coords.POINT_MODE:
             lats, lons = get_lat_lon_grids(cube)
         else:
@@ -154,7 +160,7 @@ def get_lat_lon_grids(cube):
     """
     # get the lat and lon coords (might have "grid_" at the start of the name, if rotated).
     lat_coord, lon_coord = _get_lat_lon_coords(cube)
-    cs = cube.coord_system('LatLonCS')
+    cs = cube.coord_system('CoordSystem')
 
     if lon_coord.units != 'degrees':
         lon_coord = lon_coord.unit_converted('degrees')
@@ -168,8 +174,8 @@ def get_lat_lon_grids(cube):
     lons, lats = numpy.meshgrid(lons, lats)
 
     # if the pole was rotated, then un-rotate it
-    if cs.has_rotated_pole():
-        lons, lats = unrotate_pole(lons, lats, cs.n_pole.longitude, cs.n_pole.latitude)
+    if isinstance(cs, iris.coord_systems.RotatedGeogCS):
+        lons, lats = unrotate_pole(lons, lats, cs.grid_north_pole_longitude, cs.grid_north_pole_latitude)
 
     return (lats, lons)
 
@@ -186,7 +192,7 @@ def get_lat_lon_contiguous_bounded_grids(cube):
     """
     # get the lat and lon coords (might have "grid_" at the start of the name, if rotated).
     lat_coord, lon_coord = _get_lat_lon_coords(cube)
-    cs = cube.coord_system('LatLonCS')
+    cs = cube.coord_system('CoordSystem')
 
     if lon_coord.units != 'degrees':
         lon_coord = lon_coord.unit_converted('degrees')
@@ -196,8 +202,8 @@ def get_lat_lon_contiguous_bounded_grids(cube):
     lons = lon_coord.contiguous_bounds()
     lats = lat_coord.contiguous_bounds()
     lons, lats = numpy.meshgrid(lons, lats)
-    if cs.has_rotated_pole():
-        lons, lats = iris.analysis.cartography.unrotate_pole(lons, lats, cs.n_pole.longitude, cs.n_pole.latitude)
+    if isinstance(cs, iris.coord_systems.RotatedGeogCS):
+        lons, lats = iris.analysis.cartography.unrotate_pole(lons, lats, cs.grid_north_pole_longitude, cs.grid_north_pole_latitude)
     return (lats, lons)
 
 
@@ -254,9 +260,15 @@ def area_weights(cube):
 
     """
     # Get the radius of the earth
-    latlon_cs = cube.coord_system(iris.coord_systems.LatLonCS)
-    if latlon_cs and latlon_cs.datum.is_spherical():
-        radius_of_earth = latlon_cs.datum.semi_major_axis
+    cs = cube.coord_system("CoordSystem")
+    if isinstance(cs, iris.coord_systems.GeogCS):
+        if cs.inverse_flattening != 0.0:
+            warnings.warn("Assuming spherical earth from ellipsoid.")
+        radius_of_earth = cs.semi_major_axis
+    elif isinstance(cs, iris.coord_systems.RotatedGeogCS) and cs.ellipsoid is not None:
+        if cs.ellipsoid.inverse_flattening != 0.0:
+            warnings.warn("Assuming spherical earth from ellipsoid.")
+        radius_of_earth = cs.ellipsoid.semi_major_axis
     else:
         warnings.warn("Using DEFAULT_SPHERICAL_EARTH_RADIUS.")
         radius_of_earth = DEFAULT_SPHERICAL_EARTH_RADIUS

lib/iris/analysis/interpolate.py

@@ -358,7 +358,7 @@ def regrid(source_cube, grid_cube, mode='bilinear', **kwargs):
     by the grid_cube.
 
     Fundamental input requirements:
-        1) Both cubes must have a HorizontalCS.
+        1) Both cubes must have a CoordSystem.
         2) The source 'x' and 'y' coordinates must not share data dimensions with any other coordinates.
 
     In addition, the algorithm currently used requires:
@@ -386,12 +386,12 @@ def regrid(source_cube, grid_cube, mode='bilinear', **kwargs):
 
     """
     # Condition 1
-    source_cs = source_cube.coord_system(iris.coord_systems.HorizontalCS)
-    grid_cs = grid_cube.coord_system(iris.coord_systems.HorizontalCS)
-    if source_cs is None or grid_cs is None:
-        raise ValueError("The source and grid cubes must contain a HorizontalCS.")
+    source_cs = source_cube.coord_system(iris.coord_systems.CoordSystem)
+    grid_cs = grid_cube.coord_system(iris.coord_systems.CoordSystem)
+    if (source_cs is None) != (grid_cs is None):
+        raise ValueError("The source and grid cubes must both have a CoordSystem or both have None.")
 
-    # Condition 2: We can only have one x coordinate and one y coordinate with the source HorizontalCS, and those coordinates 
+    # Condition 2: We can only have one x coordinate and one y coordinate with the source CoordSystem, and those coordinates 
     # must be the only ones occupying their respective dimension 
     source_x = source_cube.coord(axis='x', coord_system=source_cs)
     source_y = source_cube.coord(axis='y', coord_system=source_cs)