Pycarto++

CHANGES

@@ -1,3 +1,17 @@
+Release 1.0
+===========
+
+Incompatible changes
+--------------------
+* With the change to cartopy, the function iris.plot.gcm() is obsolete
+  and has been removed.
+* Three functions have been removed from iris.analysis.cartography:
+  lat_lon_range(), get_lat_lon_grids(), and
+  get_lat_lon_contiguous_grids(). They have been replaced with
+  generalised versions: xy_range(), get_xy_grids(),
+  and get_xy_contiguous_bounded_grids().
+
+
 Release 0.9 (14 Aug, 2012)
 ==========================
 

docs/iris/example_code/graphics/COP_maps.py

@@ -86,17 +86,15 @@ def main():
         plt.subplot(121)
         plt.title('HadGEM2 E1 Scenario',  fontsize=10)
         iplt.contourf(delta_e1, levels, colors=colors, linewidth=0, extend='both')
-        current_map = iplt.gcm()
-        current_map.drawcoastlines()
+        plt.gca().coastlines()
         # get the current axes' subplot for use later on
         plt1_ax = plt.gca()
 
         # Add the second subplot showing the A1B scenario
         plt.subplot(122)
         plt.title('HadGEM2 A1B-Image Scenario',  fontsize=10)
         contour_result = iplt.contourf(delta_a1b, levels, colors=colors, linewidth=0, extend='both')
-        current_map = iplt.gcm()
-        current_map.drawcoastlines()
+        plt.gca().coastlines()
         # get the current axes' subplot for use later on
         plt2_ax = plt.gca()
 

docs/iris/example_code/graphics/TEC.py

@@ -41,8 +41,8 @@ def main():
     plt.title('Total Electron Content')
     plt.xlabel('longitude / degrees')
     plt.ylabel('latitude / degrees')
-    iplt.gcm().bluemarble(zorder=-1)
-    iplt.gcm().drawcoastlines()
+    plt.gca().bluemarble()
+    plt.gca().coastlines()
     plt.show()
 
 

docs/iris/example_code/graphics/custom_file_loading.py

@@ -215,15 +215,11 @@ def main():
     # Callback shown as None to illustrate where a cube-level callback function would be used if required
     cube = iris.load_strict(fname, boundary_volc_ash_constraint, callback=None)
 
-    map = iplt.map_setup(lon_range=[-70, 20], lat_range=[20, 75], resolution='i')
-
-    map.drawcoastlines()
-
-    iplt.contourf(cube, 
-                        levels=(0.0002, 0.002, 0.004, 1),
+    iplt.map_setup(xlim=(-70, 20), ylim=(20, 75))
+    plt.gca().coastlines()
+    iplt.contourf(cube, levels=(0.0002, 0.002, 0.004, 1),
                         colors=('#80ffff', '#939598', '#e00404'),
-                        extend='max'
-                  )
+                        extend='max')
 
     time = cube.coord('time')
     time_date = time.units.num2date(time.points[0]).strftime(UTC_format)

docs/iris/example_code/graphics/global_map.py

@@ -18,7 +18,7 @@ def main():
     temperature = iris.load_strict(fname)
 
     qplt.contourf(temperature, 15)
-    iplt.gcm().drawcoastlines()
+    plt.gca().coastlines()
     plt.show()
 
 

docs/iris/example_code/graphics/lagged_ensemble.py

@@ -76,8 +76,7 @@ def main():
         cf = iplt.contourf(cube, contour_levels)
 
         # add coastlines
-        m = iplt.gcm()
-        m.drawcoastlines()
+        plt.gca().coastlines()
 
     # make an axes to put the shared colorbar in
     colorbar_axes = plt.gcf().add_axes([0.35, 0.1, 0.3, 0.05])

docs/iris/example_code/graphics/rotated_pole_mapping.py

@@ -23,36 +23,27 @@ def main():
     fname = iris.sample_data_path('rotated_pole.nc')
     temperature = iris.load_strict(fname)
 
-    # Calculate the lat lon range and buffer it by 10 degrees
-    lat_range, lon_range = iris.analysis.cartography.lat_lon_range(temperature)
-    lat_range = lat_range[0] - 10, lat_range[1] + 10
-    lon_range = lon_range[0] - 10, lon_range[1] + 10
-
-
     # Plot #1: Point plot showing data values & a colorbar
     plt.figure()
-    iplt.map_setup(temperature, lat_range=lat_range, lon_range=lon_range)
     points = qplt.points(temperature, c=temperature.data)
     cb = plt.colorbar(points, orientation='horizontal')
     cb.set_label(temperature.units)
-    iplt.gcm().drawcoastlines()
+    plt.gca().coastlines()
     plt.show()
 
 
     # Plot #2: Contourf of the point based data
     plt.figure()
-    iplt.map_setup(temperature, lat_range=lat_range, lon_range=lon_range)
     qplt.contourf(temperature, 15)
-    iplt.gcm().drawcoastlines()
+    plt.gca().coastlines()
     plt.show()
 
 
     # Plot #3: Contourf overlayed by coloured point data
     plt.figure()
-    iplt.map_setup(temperature, lat_range=lat_range, lon_range=lon_range)
     qplt.contourf(temperature)
     iplt.points(temperature, c=temperature.data)
-    iplt.gcm().drawcoastlines()
+    plt.gca().coastlines()
     plt.show()
 
 
@@ -64,10 +55,9 @@ def main():
 
     # Plot #4: Block plot
     plt.figure()
-    iplt.map_setup(temperature, lat_range=lat_range, lon_range=lon_range)
     iplt.pcolormesh(temperature)
-    iplt.gcm().bluemarble()
-    iplt.gcm().drawcoastlines()
+    plt.gca().bluemarble()
+    plt.gca().coastlines()
     plt.show()
 
 

docs/iris/src/conf.py

@@ -146,7 +146,6 @@
    'numpy': ('http://docs.scipy.org/doc/numpy/', None), 
    'scipy': ('http://docs.scipy.org/doc/scipy/reference/', None), 
    'matplotlib': ('http://matplotlib.sourceforge.net/', None),
-   'basemap': ('http://matplotlib.github.com/basemap/', None),
 } 
 
 

docs/iris/src/userguide/plotting_a_cube.rst

@@ -207,15 +207,15 @@ Plotting 2-dimensional cubes
 
 Creating maps
 -------------
-Whenever a 2D plot is created and the x and y coordinates are longitude and latitude a 
-:class:`mpl_toolkits.basemap.Basemap` instance is created which can be accessed with the :func:`iris.plot.gcm` function. 
+Whenever a 2D plot is created using an :class:`iris.coord_systems.CoordSystem` a 
+cartopy :class:`~cartopy.mpl_integration.GenericProjectionAxes` instance is created
+which can be accessed with the :func:`matplotlib.pyplot.gca` function. 
 
-Given the current map, you can draw meridians, parallels and coastlines amongst other things. 
+Given the current map, you can draw gridlines and coastlines amongst other things. 
 
 .. seealso:: 
-	:meth:`Basemap.drawmeridians() <mpl_toolkits.basemap.Basemap.drawmeridians>`, 
-	:meth:`Basemap.drawparallels() <mpl_toolkits.basemap.Basemap.drawparallels>` and 
-	:meth:`Basemap.drawcoastlines() <mpl_toolkits.basemap.Basemap.drawcoastlines>`.
+	:meth:`cartopy's gridlines() <cartopy.mpl_integration.GenericProjectionAxes.gridlines>`, 
+	:meth:`cartopy's coastlines() <cartopy.mpl_integration.GenericProjectionAxes.coastlines>`.
 
 
 Cube contour

docs/iris/src/userguide/plotting_examples/cube_blockplot.py

@@ -14,10 +14,7 @@
 # Draw the contour with 25 levels
 qplt.pcolormesh(temperature_cube)
 
-# Get the map created by pcolormesh
-current_map = iplt.gcm()
-
-# Add coastlines to the map
-current_map.drawcoastlines()
+# Add coastlines to the map created by pcolormesh
+plt.gca().coastlines()
 
 plt.show()

docs/iris/src/userguide/plotting_examples/cube_contour.py

@@ -11,11 +11,8 @@
 # Add a contour, and put the result in a variable called contour.
 contour = qplt.contour(temperature_cube)
 
-# Get the map created by contourf
-current_map = iplt.gcm()
-
-# Add coastlines to the map
-current_map.drawcoastlines()
+# Add coastlines to the map created by contour
+plt.gca().coastlines()
 
 # Add contour labels based on the contour we have just created
 plt.clabel(contour)

docs/iris/src/userguide/plotting_examples/cube_contourf.py

@@ -10,10 +10,7 @@
 # Draw the contour with 25 levels
 qplt.contourf(temperature_cube, 25)
 
-# Get the map created by contourf
-current_map = iplt.gcm()
-
-# Add coastlines to the map
-current_map.drawcoastlines()
+# Add coastlines to the map created by contourf
+plt.gca().coastlines()
 
 plt.show()

lib/iris/analysis/cartography.py

@@ -22,8 +22,9 @@
 import math
 import warnings
 
-from mpl_toolkits.basemap import pyproj
+import pyproj
 import numpy
+import cartopy.crs
 
 import iris.analysis
 import iris.coords
@@ -95,92 +96,106 @@ def _get_lat_lon_coords(cube):
     lat_coords = filter(lambda coord: "latitude" in coord.name(), cube.coords())
     lon_coords = filter(lambda coord: "longitude" in coord.name(), cube.coords())
     if len(lat_coords) > 1 or len(lon_coords) > 1:
-        raise ValueError("Calling lat_lon_range() with multiple lat or lon coords is currently disallowed")
+        raise ValueError("Calling _get_lat_lon_coords() with multiple lat or lon coords is currently disallowed")
     lat_coord = lat_coords[0]
     lon_coord = lon_coords[0]
     return (lat_coord, lon_coord)
 
 
-def lat_lon_range(cube, mode=None):
+def xy_range(cube, mode=None, projection=None):
     """
-    Return the lat & lon range of this Cube.
+    Return the x & y range of this Cube.
 
-    If the coordinate has both points & bounds, the mode keyword can be set to determine which should be
-    used in the min/max calculation. (Must be one of iris.coords.POINT_MODE or iris.coords.BOUND_MODE)
+    Args:
+
+        * cube - The cube for which to calculate xy extents.
+
+    Kwargs:
+
+        * mode - If the coordinate has bounds, use the mode keyword to specify the
+                 min/max calculation (iris.coords.POINT_MODE or iris.coords.BOUND_MODE).
+
+        * projection - Calculate the xy range in an alternative projection. 
 
     """
     # 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)
+    x_coord, y_coord = cube.coord(axis="X"), cube.coord(axis="Y")
     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 '
+    if x_coord.has_bounds() != x_coord.has_bounds():
+        raise ValueError('Cannot get the range of the x and y coordinates if they do '
                          'not have the same presence of bounds.')
 
-    if lon_coord.has_bounds():
+    if x_coord.has_bounds():
         if mode not in [iris.coords.POINT_MODE, iris.coords.BOUND_MODE]:
             raise ValueError('When the coordinate has bounds, please specify "mode".')
         _mode = mode
     else:
         _mode = iris.coords.POINT_MODE
 
+    # Get the x and y grids
     if isinstance(cs, iris.coord_systems.RotatedGeogCS):
         if _mode == iris.coords.POINT_MODE:
-            lats, lons = get_lat_lon_grids(cube)
+            x, y = get_xy_grids(cube)
         else:
-            lats, lons = get_lat_lon_contiguous_bounded_grids(cube)
+            x, y = get_xy_contiguous_bounded_grids(cube)
     else:
         if _mode == iris.coords.POINT_MODE:
-            lons = lon_coord.points
-            lats = lat_coord.points
+            x = x_coord.points
+            y = y_coord.points
         else:
-            lons = lon_coord.bounds
-            lats = lat_coord.bounds
+            x = x_coord.bounds
+            y = y_coord.bounds
+
+    if projection:
+        # source projection
+        source_cs = cube.coord_system("CoordSystem")
+        if source_cs is not None:
+            source_proj = source_cs.as_cartopy_projection()
+        else:
+            #source_proj = cartopy.crs.PlateCarree()
+            raise Exception('Unknown source coordinate system')
+
+        if source_proj != projection:
+            # TODO: Ensure there is a test for this
+            x, y = projection.transform_points(x=x, y=y, src_crs=source_proj)
 
-    if getattr(lon_coord, 'circular', False):
-        lon_range = (numpy.min(lons), numpy.min(lons) + lon_coord.units.modulus)
+    # Get the x and y range
+    if getattr(x_coord, 'circular', False):
+        x_range = (numpy.min(x), numpy.min(x) + x_coord.units.modulus)
     else: 
-        lon_range = (numpy.min(lons), numpy.max(lons))
-    
-    return ( (numpy.min(lats), numpy.max(lats)), lon_range )
+        x_range = (numpy.min(x), numpy.max(x))
+
+    y_range = (numpy.min(y), numpy.max(y))
 
+    return (x_range, y_range)
 
-def get_lat_lon_grids(cube):
+
+def get_xy_grids(cube):
     """
-    Return 2d lat and lon points in the requested coordinate system.
+    Return 2d x and y points in the native coordinate system.
     ::
 
-        lats, lons = get_lat_lon_grids(cube)
+        x, y = get_xy_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)
+    x_coord, y_coord = cube.coord(axis="X"), cube.coord(axis="Y")
     cs = cube.coord_system('CoordSystem')
 
-    if lon_coord.units != 'degrees':
-        lon_coord = lon_coord.unit_converted('degrees')
-    if lat_coord.units != 'degrees':
-        lat_coord = lat_coord.unit_converted('degrees')
-
-    lons = lon_coord.points
-    lats = lat_coord.points
+    x = x_coord.points
+    y = y_coord.points
 
     # Convert to 2 x 2d grid of data
-    lons, lats = numpy.meshgrid(lons, lats)
-
-    # if the pole was rotated, then un-rotate it
-    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)
+    x, y = numpy.meshgrid(x, y)
 
+    return (x, y)
 
-def get_lat_lon_contiguous_bounded_grids(cube):
+
+def get_xy_contiguous_bounded_grids(cube):
     """
     Return 2d lat and lon bounds.
 
@@ -190,21 +205,14 @@ def get_lat_lon_contiguous_bounded_grids(cube):
         lats, lons = cs.get_lat_lon_bounded_grids()
 
     """
-    # 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)
+    x_coord, y_coord = cube.coord(axis="X"), cube.coord(axis="Y")
     cs = cube.coord_system('CoordSystem')
 
-    if lon_coord.units != 'degrees':
-        lon_coord = lon_coord.unit_converted('degrees')
-    if lat_coord.units != 'degrees':
-        lat_coord = lat_coord.unit_converted('degrees')
-
-    lons = lon_coord.contiguous_bounds()
-    lats = lat_coord.contiguous_bounds()
-    lons, lats = numpy.meshgrid(lons, lats)
-    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)
+    x = x_coord.contiguous_bounds()
+    y = y_coord.contiguous_bounds()
+    x, y = numpy.meshgrid(x, y)
+
+    return (x, y)
 
 
 def _quadrant_area(radian_colat_bounds, radian_lon_bounds, radius_of_earth):