diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 286322697..ae32baba0 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -41,6 +41,7 @@
 import clawpack.geoclaw.data
 import six
 from six.moves import range
+from clawpack.geoclaw.data import Rearth, DEG2RAD, RAD2DEG
 
 # ==============================================================================
 #  Topography Related Functions
@@ -1672,7 +1673,383 @@ def make_shoreline_xy(self, sea_level=0):
             shoreline_xy = numpy.vstack((shoreline_xy, \
                            numpy.array([numpy.nan,numpy.nan]), c.allsegs[0][k]))
         return shoreline_xy
+       
+    
+def intersection(rect1, rect2):
+    r"""
+    Whether two domains, rect1 and rect2 intersect. If they intersect,
+    find the intersection's boundary and area.
+    :Input:
+     - *rect1* (ndarray(:))
+     - *rect2* (ndarray(:))
+    :Output:
+     - *mark* (bool) If rect1 and rect2 intersect, mark = True. Otherwise,
+       mark = False.
+     - *area* (float) The area of the intersection of rect1 and rect2.
+     - *bound* (ndarray(:)) The boundary of the intersection of rect1 and rect2.
+    """
+
+    # Set boundary for rect1
+    x1_low = rect1[0]; x1hi = rect1[1]
+    y1_low = rect1[2]; y1hi = rect1[3]
+
+    # Set boundary for rect2
+    x2low = rect2[0]; x2hi = rect2[1]
+    y2low = rect2[2]; y2hi = rect2[3]
+
+    # Boundary of the intersection part
+    xintlow = max(x1_low, x2low)
+    xinthi = min(x1hi, x2hi)
+    yintlow = max(y1_low, y2low)
+    yinthi = min(y1hi, y2hi)
+
+    # Whether rect1 and rect2 intersect
+    if xinthi > xintlow and yinthi > yintlow:
+        area = (xinthi - xintlow) * (yinthi - yintlow)
+        mark = True
+        bound = [xintlow, xinthi, yintlow, yinthi]
+    else:
+        area = 0.0
+        mark = False
+        bound = [-1, -1, -1, -1]
+
+    return mark, area, bound
+
+    
+def bilinearintegral(domain, topodomain, corners):
+    r"""
+    bilinearintegral integrates a surface over a rectangular region which is the
+    intersection of the surface defined by "domain" and the surface defined by
+    "topodomain". The surface integrated is defined by the function obtained by
+    bilinear interpolation through corners of the Cartesian grid.
+    
+    Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
+    f(xi, eta) = a * (xi - x1) / (x2 - x1) + b * (eta - y1) / (y2 - y1)
+                 + c * (xi - x1) * (eta - y1) / {(x2 - x1) * (y2 - y1)} + d,
+    to represent topo of topomain by bilinear interpolation with four corner
+    points. Then integrate f(xi, eta) on the rectangular domain which is the
+    intersection of domain and topodomain.
+    
+    :Input:
+     - *domain* (ndarray(:))
+     - *topodomain* (ndarray(:))
+     - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+    :Output:
+     - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
+       the intersection of domain and topodomain.
+    """
+
+    # Set boundary for integral
+    bound = intersection(domain, topodomain)[2]
+
+    # Find terms for the integration
+    deltax = topodomain[1] - topodomain[0]
+    deltay = topodomain[3] - topodomain[2]
+    area = (bound[3] - bound[2]) * (bound[1] - bound[0])
+    sumxi = (bound[1] + bound[0] - 2.0 * topodomain[0]) / deltax
+    sumeta = (bound[3] + bound[2] - 2.0 * topodomain[2]) / deltay
+
+    # Find coefficients of the function, f(xi, eta)
+    a = corners[1][0] - corners[0][0]
+    b = corners[0][1] - corners[0][0]
+    c = corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]
+    d = corners[0][0]
+
+    bilinearintegral = (0.5 * (a * sumxi + b * sumeta) +
+                        0.25 * c * sumxi * sumeta + d) * area
+
+    return bilinearintegral
+    
+    
+def bilinearintegral_s(domain, topodomain, corners):
+    r"""
+    bilinearintegral_s integrates a surface over a sphere which is the intersection
+    of the surface defined by "domain" and the surface defined by "topodomain".
+    The surface integrated is defined by the function obtained by bilinear
+    interpolation through corners of the polar grid.
+    
+    Suppose corners are (x1, y1), (x1, y2), (x2, y1), (x1, y2). Get a function,
+    f(theta, phi) = a * (theta - x1) + b * (phi - y1)
+                    + c * (theta - x1) * (phi - y1) + d,
+    where theta is longitute and phi is latitude, to represent the topo of
+    topomain by bilinear interpolation with four corner points. Then integrate
+    f(theta, phi) on the domain which is the intersection of the domain
+    and the topodomain on the sphere.
+    
+    :Input:
+     - *domain* (ndarray(:))
+     - *topodomain* (ndarray(:))
+     - *corners* (ndarray(:, :)) Four corner points of the topodomain.
+    :Output:
+     - *bilinearintegral* (float) The integral of the function, f(xi, eta), on
+       the intersection of domain and topodomain.
+    """
+
+    #Set boundary parameter
+    bound = intersection(domain, topodomain)[2]
+
+    # Find terms for the integration
+    xdiffhi = bound[1] - topodomain[0]
+    xdifflow = bound[0] - topodomain[0]
+    ydiffhi = bound[3] - topodomain[2]
+    ydifflow = bound[2] - topodomain[2]
+    xdiff2 = 0.5 * (xdiffhi**2 - xdifflow**2)
+    intdx = bound[1] - bound[0]
+    deltax = topodomain[1] - topodomain[0]
+    deltay = topodomain[3] - topodomain[2]
+    
+    d2r = DEG2RAD; r2d = RAD2DEG
+
+    cbsinint = (r2d * numpy.cos(d2r * bound[3]) + ydiffhi * numpy.sin(d2r * bound[3]))
+    - (r2d * numpy.cos(d2r * bound[2]) + ydifflow * numpy.sin(d2r * bound[2]))
+
+    adsinint = r2d * (numpy.sin(d2r * bound[3]) - numpy.sin(d2r * bound[2]))
+
+    # Find coefficients of the function, f(theta, phi)
+    a = (corners[1][0] - corners[0][0]) / deltax
+    b = (corners[0][1] - corners[0][0]) / deltay
+    c = (corners[1][1] - corners[1][0] - corners[0][1] + corners[0][0]) / (deltax * deltay)
+    d = corners[0][0]
+
+    bilinearintegral_s = ((a * xdiff2 + d * intdx) * adsinint +
+                          r2d * (c * xdiff2 + b * intdx) * cbsinint) * (Rearth * d2r)**2
+
+    return bilinearintegral_s
+
+
+def topointegral(domain, topoparam, z, intmethod, coordinate_system):
+    r"""
+    topointegral integrates a surface over a rectangular region which is
+    the intersection of the surface defined by "domain" and the surface
+    defined by "topoparam". The surface integrated is defined by the function
+    obtained by bilinear interpolation through nodes of the Cartesian grid.
+    
+    :Input:
+     - *domain* (ndarray(:))
+     - *topoparam* (ndarray(:)) topoparam includes topo's x and y coordinates'
+       minimum and maximum, and topo's grid size, dx and dy.
+     - *z* (ndarray(:, :)) Topo data of domain represented by topoparam.
+     - *intmethod* (int) The method of integral.
+     - *coordinate_system* (int) The type of coordinate system.
+    :Output:
+     - *theintegral* (float) The integral of a surface over a rectangular
+       region which is the intersection of the domain and "topoparam".
+    """
 
+    theintegral = 0.0
+
+    # Set topo's parameter
+    x1 = topoparam[0]; x2 = topoparam[1]
+    y1 = topoparam[2]; y2 = topoparam[3]
+    topodx = topoparam[4]; topody = topoparam[5]
+    topomx = numpy.array(z).shape[1]
+    topomy = numpy.array(z).shape[0]
+
+    # Find the intersection of the domain and topo
+    bound = intersection(domain, topoparam)[2]
+    xlow = bound[0]; xhi = bound[1]
+    ylow = bound[2]; yhi = bound[3]
+
+    if intmethod == 1:
+
+        # Find topo's start and end points for integral
+        # The x coodinate of the topo's start point
+        istart = -1
+        for i in range(topomx):
+            if (x1 + i * topodx) > xlow:
+                istart = max(i-1, 0)
+                break
+
+        # The y coodinate of the topo's start point
+        jstart = -1
+        for j in range(topomy):
+            if (y1 + j * topody)  > ylow:
+                jstart = max(j-1, 0)
+                break
+ 
+        # The x coodinate of the topo's end point
+        iend = topomx
+        for m in range(topomx):
+            if (x1 + m * topodx) >= xhi:
+                iend = m
+                break
+   
+        # The y coodinate of the topo's end point
+        jend = topomy
+        for n in range(topomy):
+            if (y1 + n * topody) >= yhi:
+                jend = n
+                break
+
+        # Prevent overflow
+        jstart = max(jstart, 0)
+        istart = max(istart, 0)
+        jend = min(jend, topomy - 1)
+        iend = min(iend, topomx - 1)
+
+        # Topo integral
+        for jj in range(jstart, jend):
+            yint1 = y1 + jj * topody
+            yint2 = y1 + (jj + 1) * topody
+
+            for ii in range(istart, iend):
+                xint1 = x1 + ii * topodx
+                xint2 = x1 + (ii + 1) * topodx
+
+                # Four corners of topo's grid
+                corners = numpy.empty((2, 2))
+                corners[0][0] = z[jj][ii]
+                corners[0][1] = z[jj+1][ii]
+                corners[1][0] = z[jj][ii+1]
+                corners[1][1] = z[jj+1][ii+1]
+
+                # Parameters of topo's grid integral
+                topointparam = [xint1, xint2, yint1, yint2]
+
+                if coordinate_system == 1:
+                    theintegral += bilinearintegral(domain, topointparam, corners)
+                elif coordinate_system == 2:
+                    theintegral += bilinearintegral_s(domain, topointparam, corners)
+                else:
+                    print("TOPOINTEGRAL: coordinate_system error")
+
+    else:
+        print("TOPOINTEGRAL: only intmethod = 1,2 is supported")
+
+    return theintegral
+
+
+def recurintegral(domain, mtopoorder, mtopofiles, m, topo):
+    r"""
+    Compute the integral of topo over the rectangle domain defined by "domain".
+    Find the index of the topo whose resolution is m by "mtopoorder". Use this
+    index to find corresponding topo object in the list "topo".
+    The main call to this recursive function has corners of a grid cell for the
+    rectangle and m = 1 in order to compute the integral over the cell using all
+    topo objects.
+    The recursive strategy is to first compute the integral using only the topo
+    object of m+1 resolution. Then apply corrections due to adding m resolution
+    topo object.
+     
+    Corrections are needed if the new topo object intersects the grid cell.
+    Let the intersection be represented by mark2[2]. Two corrections are needed.
+    First, subtract out the integral over the rectangle "mark2[2]" computed with
+    topo objects mtopoorder(mtopofiles) to mtopoorder(m+1), and then adding in
+    the integral over this same region using the topo object mtopoorder(m).
+    
+    :Input:
+     - *domain* (ndarray(:)) The specific surface where the integral is computed.
+     - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+     - *mtopofiles* (int) The number of the topo objects.
+     - *m* (int) Resolution of the topo objects.
+     - *topo* (list) The list of topo objects.
+    :Output:
+     - *integral* (float) Integral.
+    """
+    
+    # The index of the topo object whose resolution is m
+    mfid = mtopoorder[m]
+    
+    # The parameters of the m resolution topo object
+    topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
+                      topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
+    
+    # Innermost step of recursion reaches this point, only using coarsest topo grid
+    # compute directly
+    if m == mtopofiles - 1:
+        mark1 = intersection(domain, topoparam_mfid)
+        if mark1[0]:
+            integral = topointegral(mark1[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+        else:
+            integral = 0.0
+    else:
+        int1 = recurintegral(domain, mtopoorder, mtopofiles, m+1, topo)
+
+        # Whether new topo object intersects the grid cell
+        mark2 = intersection(domain, topoparam_mfid)
+        
+        # The new topo object intersects the grid cell. Corrections are needed
+        if mark2[1] > 0:
+            int2 = recurintegral(mark2[2], mtopoorder, mtopofiles, m+1, topo)
+            int3 = topointegral(mark2[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+            integral = int1 - int2 + int3
+        else:
+            integral = int1
+    return integral
+        
+
+def cellintegral(cell, mtopoorder, mtopofiles, topo):
+    r"""
+    Compute the integral on the surface defined by "cell" with the topo objects
+    list "topo".
+    
+    :Input:
+    - *cell* (ndarray(:)) The boundary of specific surface where the integral
+      is computed.
+    - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+    - *mtopofiles* (int) The number of the topo objects.
+    - *topo* (list) The list of topo objects.
+    :Output:
+    - *topoint* (float) The integral on the cell, "cell".
+    """
+    
+    # Initialization
+    topoint = 0.0
+
+    for m in range(mtopofiles):
+    
+        # The index of the topo file object resolution is m
+        mfid = mtopoorder[m]
+        cellarea = (cell[1] - cell[0]) * (cell[3] - cell[2])
+        
+        # The parameters of this topo object
+        topoparam_mfid = [topo[mfid].extent[0], topo[mfid].extent[1], topo[mfid].extent[2],
+                          topo[mfid].extent[3], topo[mfid].delta[0], topo[mfid].delta[1]]
+        
+        # Whether m-th topo interects with cell
+        mark = intersection(cell, topoparam_mfid)
+        if mark[0]:
+            
+            # Whether the cell is completely overlapped by m-th topo
+            if mark[1] == cellarea:
+                topoint = topoint + topointegral(mark[2], topoparam_mfid, topo[mfid].Z, 1, 1)
+                return topoint
+            else:
+                
+                # If the cell is not completely overlapped by m-th topo, use recursion
+                # to compute the integral on the cell
+                topoint = recurintegral(cell, mtopoorder, mtopofiles, m, topo)
+                return topoint
+    return topoint
+
+
+def cell_average_patch(patch, mtopoorder, mtopofiles, topo):
+    r"""
+    Compute every cell's average of the specific patch given by "patch".
+    :Input:
+    - *patch* (Dictionary) Patch.
+    - *mtopoorder* (ndarray(:)) The order of the topo objects' resolutions.
+    - *mtopofiles* (int) The number of the topo objects.
+    - *topo* (list) The list of topo objects.
+    :Output:
+    - *cell_average* (ndarray(:, :)) Patch's every cell average.
+    """
+
+    x_num = len(patch['x']) - 1
+    y_num = len(patch['y']) - 1
+    cell_average = numpy.empty((y_num, x_num))
+
+    for i in range(y_num):
+        for j in range(x_num):
+            
+            # Set cell's parameter
+            cell = [patch['x'][j], patch['x'][j+1], patch['y'][i], patch['y'][i+1]]
+            area = (patch['x'][j+1] - patch['x'][j]) * (patch['y'][i+1] - patch['y'][i])
+            
+            # Calculate every cell average
+            cell_average[i, j] = cellintegral(cell, mtopoorder, mtopofiles, topo) / area
+            
+    return cell_average
 
 
 # Define convenience dictionary of URLs for some online DEMs in netCDF form:
diff --git a/tests/test_topotools.py b/tests/test_topotools.py
index ed07f583f..f043cad64 100644
--- a/tests/test_topotools.py
+++ b/tests/test_topotools.py
@@ -20,7 +20,10 @@
 
 import clawpack.geoclaw.topotools as topotools
 import clawpack.clawutil.data
+import numpy.testing as npt
 from six.moves import range
+from clawpack.geoclaw.topotools import Topography
+from scipy.integrate import dblquad, nquad
 
 # Set local test directory to get local files
 testdir = os.path.dirname(__file__)
@@ -427,6 +430,122 @@ def plot_kahului():
     plt.savefig(fname)
     print("Created ",fname)
 
+def integral_topotool(func, mfile, funcflag):
+    r"""
+    integral_topotool is used for testing a set of functions which compute
+    cell integrals and cell averages in topotools.py.
+    
+    :Input:
+     - *func* (function).
+     - *mfile* (int) The number of the topo objects.
+     - *funcflag* (bool) Whether the function, "func" is discontinuous.
+    """
+
+    dx = []
+    data = []
+    topo = []
+    mtopoorder = []
+    
+    # Set boundary for the coarsest topo
+    xlow = 5; xhi = 8
+    xarray = numpy.linspace(xlow, xhi, mfile * 5)
+    ylow = 2; yhi = 5
+    yarray = numpy.linspace(ylow, yhi, mfile * 5)
+    
+    # Set a set of topo data
+    for i in range(mfile):
+        if i == 0:
+            
+            # Set the topo data set which covers the whole patch
+            x = xarray; y = yarray
+            dx.append(x[1] - x[0])
+        else:
+            
+            # Set x coordinate of the topo data set
+            n1 = int(3 * i)
+            n2 = int(3 * i + 1.2 * mfile)
+            x1 = xarray[n1]; x2 = xarray[n2]
+            dx.append(dx[i-1] - 0.0012)
+            mx = int((x2 - x1) / dx[i])
+            x2 = x1 + dx[i] * (mx - 1)
+            x = numpy.linspace(x1, x2, mx)
+
+            # Set y coordinate of the topo data set
+            m1 = int(3 * i)
+            m2 = int(3 * i + 1.2 * mfile)
+            y1 = yarray[m1]; y2 = yarray[m2]
+            my = int((y2 - y1) / dx[i])
+            y2 = y1 + dx[i] * (my - 1)
+            y = numpy.linspace(y1, y2, my)
+        
+        # Set Topography objects parameters
+        topo1 = Topography()
+        
+        # Whether the function is discontinuous function
+        if funcflag:
+            z = numpy.empty((len(x), len(y)))
+            for m in range(len(x)):
+                for n in range(len(y)):
+                    z[m][n] = func(x[m], y[n])
+            topo1.Z = z
+        else:
+            topo_x, topo_y = numpy.meshgrid(x, y)
+            topo1.Z = func(topo_x, topo_y)
+        topo1.x = x; topo1.y = y
+        mtopoorder.append(mfile - 1 - i)
+        topo.append(topo1)
+    
+    # Set patch data
+    patch_x = numpy.linspace(xlow + 0.1, xhi - 0.1, 5)
+    patch_y = numpy.linspace(ylow + 0.1, yhi - 0.1, 4)
+    patch = {'x':patch_x, 'y':patch_y}
+    
+    # Accurate cell value
+    real_value = numpy.empty((len(patch['y']) - 1, len(patch['x']) - 1))
+    for i in range(len(patch['y']) - 1):
+        for j in range(len(patch['x']) - 1):
+            area = (patch['x'][j+1] - patch['x'][j]) * (patch['y'][i+1] - patch['y'][i])
+            if funcflag:
+                options = {'limit':1000}
+                real_value[i][j] = nquad(func=func, ranges=[[patch['x'][j], patch['x'][j+1]],
+                                         [patch['y'][i], patch['y'][i+1]]], args=None,
+                                         opts=[options,options])[0] / float(area)
+            else:
+                real_value[i][j] = dblquad(func=func, a=patch['y'][i], b=patch['y'][i+1],
+                                           gfun=patch['x'][j],
+                                           hfun=patch['x'][j+1])[0] / float(area)
+        
+    
+    # Cell value calculated by functions
+    calculated_value = topotools.cell_average_patch(patch, mtopoorder, mfile, topo)
+        
+    # Whether calculated value achieve the expected resolution
+    npt.assert_almost_equal(real_value, calculated_value, decimal=3)
+    return None
+
+
+def testcontinuous():
+    r"""
+    testcontinuous is used to test continuous functions with `integral_topotool`.
+    """
+    
+    integral_topotool(lambda x, y: numpy.sin(x**2 + x * y) + x + y, 20, False)
+    integral_topotool(lambda x, y: numpy.exp(0.5 * x + 0.1 * y), 20, False)
+    integral_topotool(lambda x, y: x**2 + x * y + y**2 + 1, 20, False)
+    return None
+    
+    
+def testdiscontinuous():
+    r"""
+    testdiscontinuous is used to test discontinuous function with `integral_topotool`.
+    """
+    
+    def fun(x, y):
+        for i in range(20):
+            if int(x * 100) % 20 == i:
+                return 9.99 + 0.001 * i
+    integral_topotool(fun, 20, True)
+    return None
 
 if __name__ == "__main__":
     if len(sys.argv) > 1: