def density_grid(*args, **kwargs):
    """
    Estimates point density of the given linear orientation measurements
    (Interpreted as poles, lines, rakes, or "raw" longitudes and latitudes
    based on the `measurement` keyword argument.). Returns a regular (in
    lat-long space) grid of density estimates over a hemispherical surface.

    Parameters
    ----------

    *args : 2 or 3 sequences of measurements
        By default, this will be expected to be ``strike`` & ``dip``, both
        array-like sequences representing poles to planes.  (Rake measurements
        require three parameters, thus the variable number of arguments.) The
        ``measurement`` kwarg controls how these arguments are interpreted.

    measurement : string, optional
        Controls how the input arguments are interpreted. Defaults to
        ``"poles"``.  May be one of the following:

            ``"poles"`` : strikes, dips
                Arguments are assumed to be sequences of strikes and dips of
                planes. Poles to these planes are used for contouring.
            ``"lines"`` : plunges, bearings
                Arguments are assumed to be sequences of plunges and bearings
                of linear features.
            ``"rakes"`` : strikes, dips, rakes
                Arguments are assumed to be sequences of strikes, dips, and
                rakes along the plane.
            ``"radians"`` : lon, lat
                Arguments are assumed to be "raw" longitudes and latitudes in
                the stereonet's underlying coordinate system.

    method : string, optional
        The method of density estimation to use. Defaults to
        ``"exponential_kamb"``. May be one of the following:

        ``"exponential_kamb"`` : Kamb with exponential smoothing
            A modified Kamb method using exponential smoothing [1]_. Units are
            in numbers of standard deviations by which the density estimate
            differs from uniform.
        ``"linear_kamb"`` : Kamb with linear smoothing
            A modified Kamb method using linear smoothing [1]_.  Units are in
            numbers of standard deviations by which the density estimate
            differs from uniform.
        ``"kamb"`` : Kamb with no smoothing
            Kamb's method [2]_ with no smoothing. Units are in numbers of
            standard deviations by which the density estimate differs from
            uniform.
        ``"schmidt"`` : 1% counts
            The traditional "Schmidt" (a.k.a. 1%) method. Counts points within
            a counting circle comprising 1% of the total area of the
            hemisphere. Does not take into account sample size.  Units are in
            points per 1% area.

    sigma : int or float, optional
        The number of standard deviations defining the expected number of
        standard deviations by which a random sample from a uniform
        distribution of points would be expected to vary from being evenly
        distributed across the hemisphere.  This controls the size of the
        counting circle, and therefore the degree of smoothing.  Higher sigmas
        will lead to more smoothing of the resulting density distribution. This
        parameter only applies to Kamb-based methods.  Defaults to 3.

    gridsize : int or 2-item tuple of ints, optional
        The size of the grid that the density is estimated on. If a single int
        is given, it is interpreted as an NxN grid. If a tuple of ints is given
        it is interpreted as (nrows, ncols).  Defaults to 100.

    weights : array-like, optional
        The relative weight to be applied to each input measurement. The array
        will be normalized to sum to 1, so absolute value of the weights do not
        affect the result. Defaults to None.

    Returns
    -------
    xi, yi, zi : 2D arrays
        The longitude, latitude and density values of the regularly gridded
        density estimates. Longitude and latitude are in radians.

    See Also
    ---------
    mplstereonet.StereonetAxes.density_contourf
    mplstereonet.StereonetAxes.density_contour

    References
    ----------
    .. [1] Vollmer, 1995. C Program for Automatic Contouring of Spherical
       Orientation Data Using a Modified Kamb Method. Computers &
       Geosciences, Vol. 21, No. 1, pp. 31--49.

    .. [2] Kamb, 1959. Ice Petrofabric Observations from Blue Glacier,
       Washington, in Relation to Theory and Experiment. Journal of
       Geophysical Research, Vol. 64, No. 11, pp. 1891--1909.
    """
    def do_nothing(x, y):
        return x, y
    measurement = kwargs.get('measurement', 'poles')
    gridsize = kwargs.get('gridsize', 100)
    weights = kwargs.get('weights', None)
    try:
        gridsize = int(gridsize)
        gridsize = (gridsize, gridsize)
    except TypeError:
        pass

    func = {'poles':stereonet_math.pole,
            'lines':stereonet_math.line,
            'rakes':stereonet_math.rake,
            'radians':do_nothing}[measurement]
    lon, lat = func(*args)

    method = kwargs.get('method', 'exponential_kamb')
    sigma = kwargs.get('sigma', 3)
    func = {'linear_kamb':_linear_inverse_kamb,
            'square_kamb':_square_inverse_kamb,
            'schmidt':_schmidt_count,
            'kamb':_kamb_count,
            'exponential_kamb':_exponential_kamb,
            }[method]
    lon, lat, z = _count_points(lon, lat, func, sigma, gridsize, weights)
    return lon, lat, z
Estimates point density of the given linear orientation measurements
    (Interpreted as poles, lines, rakes, or "raw" longitudes and latitudes
    based on the `measurement` keyword argument.). Returns a regular (in
    lat-long space) grid of density estimates over a hemispherical surface.

    Parameters
    ----------

    *args : 2 or 3 sequences of measurements
        By default, this will be expected to be ``strike`` & ``dip``, both
        array-like sequences representing poles to planes.  (Rake measurements
        require three parameters, thus the variable number of arguments.) The
        ``measurement`` kwarg controls how these arguments are interpreted.

    measurement : string, optional
        Controls how the input arguments are interpreted. Defaults to
        ``"poles"``.  May be one of the following:

            ``"poles"`` : strikes, dips
                Arguments are assumed to be sequences of strikes and dips of
                planes. Poles to these planes are used for contouring.
            ``"lines"`` : plunges, bearings
                Arguments are assumed to be sequences of plunges and bearings
                of linear features.
            ``"rakes"`` : strikes, dips, rakes
                Arguments are assumed to be sequences of strikes, dips, and
                rakes along the plane.
            ``"radians"`` : lon, lat
                Arguments are assumed to be "raw" longitudes and latitudes in
                the stereonet's underlying coordinate system.

    method : string, optional
        The method of density estimation to use. Defaults to
        ``"exponential_kamb"``. May be one of the following:

        ``"exponential_kamb"`` : Kamb with exponential smoothing
            A modified Kamb method using exponential smoothing [1]_. Units are
            in numbers of standard deviations by which the density estimate
            differs from uniform.
        ``"linear_kamb"`` : Kamb with linear smoothing
            A modified Kamb method using linear smoothing [1]_.  Units are in
            numbers of standard deviations by which the density estimate
            differs from uniform.
        ``"kamb"`` : Kamb with no smoothing
            Kamb's method [2]_ with no smoothing. Units are in numbers of
            standard deviations by which the density estimate differs from
            uniform.
        ``"schmidt"`` : 1% counts
            The traditional "Schmidt" (a.k.a. 1%) method. Counts points within
            a counting circle comprising 1% of the total area of the
            hemisphere. Does not take into account sample size.  Units are in
            points per 1% area.

    sigma : int or float, optional
        The number of standard deviations defining the expected number of
        standard deviations by which a random sample from a uniform
        distribution of points would be expected to vary from being evenly
        distributed across the hemisphere.  This controls the size of the
        counting circle, and therefore the degree of smoothing.  Higher sigmas
        will lead to more smoothing of the resulting density distribution. This
        parameter only applies to Kamb-based methods.  Defaults to 3.

    gridsize : int or 2-item tuple of ints, optional
        The size of the grid that the density is estimated on. If a single int
        is given, it is interpreted as an NxN grid. If a tuple of ints is given
        it is interpreted as (nrows, ncols).  Defaults to 100.

    weights : array-like, optional
        The relative weight to be applied to each input measurement. The array
        will be normalized to sum to 1, so absolute value of the weights do not
        affect the result. Defaults to None.

    Returns
    -------
    xi, yi, zi : 2D arrays
        The longitude, latitude and density values of the regularly gridded
        density estimates. Longitude and latitude are in radians.

    See Also
    ---------
    mplstereonet.StereonetAxes.density_contourf
    mplstereonet.StereonetAxes.density_contour

    References
    ----------
    .. [1] Vollmer, 1995. C Program for Automatic Contouring of Spherical
       Orientation Data Using a Modified Kamb Method. Computers &
       Geosciences, Vol. 21, No. 1, pp. 31--49.

    .. [2] Kamb, 1959. Ice Petrofabric Observations from Blue Glacier,
       Washington, in Relation to Theory and Experiment. Journal of
       Geophysical Research, Vol. 64, No. 11, pp. 1891--1909.