Cell area array broken for Voronoi grids
As reported by @SiccarPoint, the `cell_areas` array is broken. When trying to use the `cell_areas` attribute for Voronoi grids the following error is raised,

``` python
AttributeError looking for _setup_cell_areas_array
```

diff --git a/landlab/components/flow_accum/flow_accumulation.py b/landlab/components/flow_accum/flow_accumulation.py
index c38b51a9f..5a4d6f874 100644
--- a/landlab/components/flow_accum/flow_accumulation.py
+++ b/landlab/components/flow_accum/flow_accumulation.py
@@ -44,8 +44,8 @@ class AccumFlow(object):
             sorted_flowdirs = (data.flowdirs[active_cell_ids])[height_order_active_cells]
         except:
             six.print_('Flow directions could not be sorted by elevation. Does the data object contain the flow direction vector?')
-        #print grid.cell_areas
-        self.flow_accum_by_area[active_cell_ids] = grid.cell_areas #This is only the active nodes == cells by definition
+        # print grid.area_of_cell
+        self.flow_accum_by_area[active_cell_ids] = grid.area_of_cell # This is only the active nodes == cells by definition
 
         #print len(height_order_active_cells), len(sorted_flowdirs), len(self.flow_accum_by_area)
         #print height_order_active_cells
diff --git a/landlab/components/flow_accum/flow_accumulation2.py b/landlab/components/flow_accum/flow_accumulation2.py
index f61164aeb..ea9c1266a 100644
--- a/landlab/components/flow_accum/flow_accumulation2.py
+++ b/landlab/components/flow_accum/flow_accumulation2.py
@@ -41,8 +41,8 @@ class AccumFlow(object):
             sorted_flowdirs = (flowdirs[active_cell_ids])[height_order_active_cells]
         except:
             print('Flow directions could not be sorted by elevation. Does the data object contain the flow direction vector?')
-        #print grid.cell_areas
-        self.flow_accum_by_area[active_cell_ids] = grid.cell_areas #This is only the active nodes == cells by definition
+        # print grid.area_of_cell
+        self.flow_accum_by_area[active_cell_ids] = grid.area_of_cell # This is only the active nodes == cells by definition
 
         #print len(height_order_active_cells), len(sorted_flowdirs), len(self.flow_accum_by_area)
         #print height_order_active_cells
diff --git a/landlab/components/stream_power/sed_flux_dep_incision.py b/landlab/components/stream_power/sed_flux_dep_incision.py
index e61fb5969..a3726fd22 100644
--- a/landlab/components/stream_power/sed_flux_dep_incision.py
+++ b/landlab/components/stream_power/sed_flux_dep_incision.py
@@ -268,8 +268,8 @@ class SedDepEroder(object):
             print("Threshold derived from grain size and Shields number is: ", self.thresh)
 
         self.cell_areas = np.empty(grid.number_of_nodes)
-        self.cell_areas.fill(np.mean(grid.cell_areas))
-        self.cell_areas[grid.node_at_cell] = grid.cell_areas
+        self.cell_areas.fill(np.mean(grid.area_of_cell))
+        self.cell_areas[grid.node_at_cell] = grid.area_of_cell
         #new 11/12/14
         self.point6onelessb = 0.6*(1.-self._b)
         self.shear_stress_prefactor = self.fluid_density*self.g*(self.mannings_n/self.k_w)**0.6
diff --git a/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional.py b/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional.py
index 33392f248..67ce964b5 100644
--- a/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional.py
+++ b/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional.py
@@ -225,8 +225,8 @@ class TransportLimitedEroder(object):
             self.shields_prefactor_noD = 1./((self.sed_sensity-self.fluid_density)*self.g)
 
         self.cell_areas = np.empty(grid.number_of_nodes)
-        self.cell_areas.fill(np.mean(grid.cell_areas))
-        self.cell_areas[grid.node_at_cell] = grid.cell_areas
+        self.cell_areas.fill(np.mean(grid.area_of_cell))
+        self.cell_areas[grid.node_at_cell] = grid.area_of_cell
 
 
     def sed_capacity_equation(self, grid, shields_stress, slopes_at_nodes=None, areas_at_node=None):
diff --git a/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_alt.py b/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_alt.py
index 20e7ee52d..b916c4de4 100644
--- a/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_alt.py
+++ b/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_alt.py
@@ -245,8 +245,8 @@ class TransportLimitedEroder(object):
         self.diffusivity_power_on_A = 0.9*self._c*(1.-self._b) #i.e., q/D**(1/6)
 
         self.cell_areas = np.empty(grid.number_of_nodes)
-        self.cell_areas.fill(np.mean(grid.cell_areas))
-        self.cell_areas[grid.node_at_cell] = grid.cell_areas
+        self.cell_areas.fill(np.mean(grid.area_of_cell))
+        self.cell_areas[grid.node_at_cell] = grid.area_of_cell
         self.dx2 = grid.dx ** 2
         self.dy2 = grid.dy ** 2
         self.bad_neighbor_mask = np.equal(grid.get_active_neighbors_at_node(bad_index=-1),-1)
diff --git a/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_v3.py b/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_v3.py
index 6851dd9cf..939e97566 100644
--- a/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_v3.py
+++ b/landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_v3.py
@@ -245,8 +245,8 @@ class TransportLimitedEroder(object):
             self.cell_areas = grid.dx * grid.dy
         else:
             self.cell_areas = np.empty(grid.number_of_nodes)
-            self.cell_areas.fill(np.mean(grid.cell_areas))
-            self.cell_areas[grid.node_at_cell] = grid.cell_areas
+            self.cell_areas.fill(np.mean(grid.area_of_cell))
+            self.cell_areas[grid.node_at_cell] = grid.area_of_cell
         self.bad_neighbor_mask = np.equal(grid.get_active_neighbors_at_node(bad_index=-1),-1)
 
         self.routing_code = """
diff --git a/landlab/components/transport_limited_fluvial/tl_fluvial_polydirectional.py b/landlab/components/transport_limited_fluvial/tl_fluvial_polydirectional.py
index cac37a920..a57a6f128 100644
--- a/landlab/components/transport_limited_fluvial/tl_fluvial_polydirectional.py
+++ b/landlab/components/transport_limited_fluvial/tl_fluvial_polydirectional.py
@@ -245,8 +245,8 @@ class TransportLimitedEroder(object):
         self.diffusivity_power_on_A = 0.9*self._c*(1.-self._b) #i.e., q/D**(1/6)
 
         self.cell_areas = np.empty(grid.number_of_nodes)
-        self.cell_areas.fill(np.mean(grid.cell_areas))
-        self.cell_areas[grid.node_at_cell] = grid.cell_areas
+        self.cell_areas.fill(np.mean(grid.area_of_cell))
+        self.cell_areas[grid.node_at_cell] = grid.area_of_cell
         self.dx2 = grid.dx ** 2
         self.dy2 = grid.dy ** 2
         self.bad_neighbor_mask = np.equal(grid.get_active_neighbors_at_node(bad_index=-1),-1)
diff --git a/landlab/grid/base.py b/landlab/grid/base.py
index 217133378..0ce657f9c 100644
--- a/landlab/grid/base.py
+++ b/landlab/grid/base.py
@@ -1588,30 +1588,6 @@ class ModelGrid(ModelDataFields):
                                                          active_link_flux,
                                                          out=out)
 
-    @property
-    @make_return_array_immutable
-    def cell_areas(self):
-        """Cell areas.
-
-        Returns
-        -------
-        ndarray
-            Array of grid-cell areas.
-
-        Notes
-        -----
-
-        Sometimes it may make sense for a grid to not always calculate
-        its cell areas but, instead, only calculate them once they are
-        required. In such cases, the grid class must implement a
-        _setup_cell_areas_array method, which will be called the first
-        time cell areas are requested.
-        """
-        try:
-            return self._cell_areas
-        except AttributeError:
-            return self._setup_cell_areas_array()
-
     @property
     @make_return_array_immutable
     def forced_cell_areas(self):
@@ -1648,14 +1624,14 @@ class ModelGrid(ModelDataFields):
         ever need to be called.
         """
         self._forced_cell_areas = numpy.empty(self.number_of_nodes)
-        mean_cell_area = numpy.mean(self.active_cell_areas)
+        mean_cell_area = numpy.mean(self.area_of_cell)
         self._forced_cell_areas.fill(mean_cell_area)
         cell_node_ids = np.where(self.status_at_node != CLOSED_BOUNDARY)[0]
         try:
-            self._forced_cell_areas[cell_node_ids] = self.cell_areas
+            self._forced_cell_areas[cell_node_ids] = self.area_of_cell
         except AttributeError:
             # in the case of the Voronoi
-            self._forced_cell_areas[cell_node_ids] = self.active_cell_areas
+            self._forced_cell_areas[cell_node_ids] = self.area_of_cell
         return self._forced_cell_areas
 
     def get_active_link_connecting_node_pair(self, node1, node2):
@@ -1687,15 +1663,20 @@ class ModelGrid(ModelDataFields):
         return numpy.array([active_link])
 
     @property
-    def active_link_length(self):
-        """Get array of lengths of active links.
+    @make_return_array_immutable
+    def area_of_cell(self):
+        """Get areas of grid cells.
 
         Returns
         -------
         ndarray
-            Lengths of active links, in ID order.
+            Areas of cells, in ID order.
         """
-        return self.link_length[self.active_link_ids]
+        # return self._area_of_cell
+        try:
+            return self._area_of_cell
+        except AttributeError:
+            return self._setup_cell_areas_array()
 
     @property
     def link_length(self):
diff --git a/landlab/grid/grid_funcs.py b/landlab/grid/grid_funcs.py
index 5947aa85b..8fc1835df 100644
--- a/landlab/grid/grid_funcs.py
+++ b/landlab/grid/grid_funcs.py
@@ -27,10 +27,10 @@ def resolve_values_on_active_links(grid, active_link_values):
     return (
         np.multiply(((grid.node_x[grid.activelink_tonode] -
                       grid.node_x[grid.activelink_fromnode]) /
-                     grid.active_link_length), active_link_values),
+                     grid.link_length[grid.active_links]), active_link_values),
         np.multiply(((grid.node_y[grid.activelink_tonode] -
                       grid.node_y[grid.activelink_fromnode]) /
-                     grid.active_link_length), active_link_values))
+                     grid.link_length[grid.active_links]), active_link_values))
 
 
 def resolve_values_on_links(grid, link_values):
@@ -137,6 +137,6 @@ def calculate_flux_divergence_at_nodes(grid, active_link_flux, out=None):
 
     # Now divide by cell areas ... where there are core cells.
     node_at_active_cell = grid.node_at_cell[grid.core_cells]
-    net_unit_flux[node_at_active_cell] /= grid.cell_areas[grid.core_cells]
+    net_unit_flux[node_at_active_cell] /= grid.area_of_cell[grid.core_cells]
 
     return net_unit_flux
diff --git a/landlab/grid/hex.py b/landlab/grid/hex.py
index 45de50568..cfeba16b2 100644
--- a/landlab/grid/hex.py
+++ b/landlab/grid/hex.py
@@ -187,9 +187,9 @@ class HexModelGrid(VoronoiDelaunayGrid):
 
             A = 3 dx^2 / 2 \\sqrt{3} \\approx 0.866 dx^2
         """
-        self._cell_areas = 0.8660254 * self._dx**2 + \\
-            numpy.zeros(self.number_of_cells)
-        return self._cell_areas
+        self._area_of_cell = (0.8660254 * self._dx ** 2 +
+                              numpy.zeros(self.number_of_cells))
+        return self._area_of_cell
 
     @staticmethod
     def make_hex_points_horizontal_hex(num_rows, base_num_cols, dxh):
diff --git a/landlab/grid/tests/test_raster_grid/test_cell_areas.py b/landlab/grid/tests/test_raster_grid/test_cell_areas.py
index f9cce1adb..b1d83a760 100644
--- a/landlab/grid/tests/test_raster_grid/test_cell_areas.py
+++ b/landlab/grid/tests/test_raster_grid/test_cell_areas.py
@@ -11,31 +11,31 @@ def setup():
 
 @with_setup(setup)
 def test_unit_grid_all_cells():
-    assert_array_equal(rmg.cell_areas, np.ones(4))
+    assert_array_equal(rmg.area_of_cell, np.ones(4))
 
 
 @with_setup(setup)
 def test_unit_grid_one_cell():
-    assert_array_equal(rmg.cell_areas[0], 1.)
+    assert_array_equal(rmg.area_of_cell[0], 1.)
 
 
 @with_setup(setup)
 def test_unit_grid_last_cell():
-    assert_array_equal(rmg.cell_areas[-1], 1.)
+    assert_array_equal(rmg.area_of_cell[-1], 1.)
 
 
 @with_setup(setup)
 @raises(IndexError)
 def test_out_of_range():
-    rmg.cell_areas[5]
+    rmg.area_of_cell[5]
 
 
 @with_setup(setup)
 @raises(ValueError)
 def test_is_immutable():
-    rmg.cell_areas[0] = 0.
+    rmg.area_of_cell[0] = 0.
 
 
 def test_all_cells_with_spacing():
     rmg = RasterModelGrid(4, 4, 10.)
-    assert_array_equal(rmg.cell_areas, 100 * np.ones(4))
+    assert_array_equal(rmg.area_of_cell, 100 * np.ones(4))
diff --git a/landlab/grid/tests/test_raster_grid/test_link_length.py b/landlab/grid/tests/test_raster_grid/test_link_length.py
index 1a3b10e19..624654cf8 100644
--- a/landlab/grid/tests/test_raster_grid/test_link_length.py
+++ b/landlab/grid/tests/test_raster_grid/test_link_length.py
@@ -37,5 +37,5 @@ def test_link_length():
 
 @with_setup(setup_grid)
 def test_active_link_length():
-    assert_array_equal(_RMG.active_link_length,
+    assert_array_equal(_RMG.link_length[_RMG.active_links],
                        [3.] * 9  + [4.] * 8 )
diff --git a/landlab/grid/voronoi.py b/landlab/grid/voronoi.py
index 876034013..c162c857f 100644
--- a/landlab/grid/voronoi.py
+++ b/landlab/grid/voronoi.py
@@ -183,11 +183,11 @@ class VoronoiDelaunayGrid(ModelGrid):
         # each active cell.
         vor = Voronoi(pts)
         self.vor = vor
-        self.active_cell_areas = numpy.zeros(self.number_of_active_cells)
+        self._area_of_cell = numpy.zeros(self.number_of_active_cells)
         for node in self._node_at_cell:
             xv = vor.vertices[vor.regions[vor.point_region[node]], 0]
             yv = vor.vertices[vor.regions[vor.point_region[node]], 1]
-            self.active_cell_areas[self.cell_at_node[node]] = (
+            self._area_of_cell[self.cell_at_node[node]] = (
                 simple_poly_area(xv, yv))
 
         # LINKS: Construct Delaunay triangulation and construct lists of link
diff --git a/landlab/plot/channel_profile.py b/landlab/plot/channel_profile.py
index 8a7d24aac..c77eecad6 100644
--- a/landlab/plot/channel_profile.py
+++ b/landlab/plot/channel_profile.py
@@ -27,7 +27,7 @@ except ImportError:
 
 def channel_nodes(grid, steepest_nodes, drainage_area, flow_receiver, number_of_channels=1, threshold=None):
     if threshold == None:
-        threshold = 2. * numpy.amin(grid.cell_areas)
+        threshold = 2. * numpy.amin(grid.area_of_cell)
     boundary_nodes = grid.boundary_nodes
     #top_two_pc = len(boundary_nodes)//50
     #starting_nodes = boundary_nodes[numpy.argsort(drainage_area[boundary_nodes])[-top_two_pc:]]
@@ -125,7 +125,7 @@ def analyze_channel_network_and_plot(grid, elevations='topographic__elevation',
         assert len(
             starting_nodes) == number_of_channels, "Length of starting_nodes must equal the number_of_channels!"
         if threshold == None:
-            threshold = 2. * numpy.amin(grid.cell_areas)
+            threshold = 2. * numpy.amin(grid.area_of_cell)
         profile_IDs = []
         for i in starting_nodes:
             j = i
['landlab/components/flow_accum/flow_accumulation.py', 'landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional.py', 'landlab/grid/voronoi.py', 'landlab/components/transport_limited_fluvial/tl_fluvial_polydirectional.py', 'landlab/grid/hex.py', 'landlab/grid/grid_funcs.py', 'landlab/grid/tests/test_raster_grid/test_link_length.py', 'landlab/components/stream_power/sed_flux_dep_incision.py', 'landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_alt.py', 'landlab/components/transport_limited_fluvial/tl_fluvial_monodirectional_v3.py', 'landlab/components/flow_accum/flow_accumulation2.py', 'landlab/grid/base.py', 'landlab/grid/tests/test_raster_grid/test_cell_areas.py', 'landlab/plot/channel_profile.py']