How to use the wradlib.georef.projection.get_earth_radius function in wradlib

xyz : :class:`numpy:numpy.ndarray`
        Array of shape (..., 3). Contains cartesian coordinates.
    rad : osr object
        Destination Spatial Reference System (Projection).
        Defaults to wgs84 (epsg 4326).
    """
    # if site altitude is present, use it, else assume it to be zero
    try:
        centalt = sitecoords[2]
    except IndexError:
        centalt = 0.

    # if no radius is given, get the approximate radius of the WGS84
    # ellipsoid for the site's latitude
    if re is None:
        re = projection.get_earth_radius(sitecoords[1])
        # Set up aeqd-projection sitecoord-centered, wgs84 datum and ellipsoid
        # use world azimuthal equidistant projection
        projstr = ('+proj=aeqd +lon_0={lon:f} +x_0=0 +y_0=0 +lat_0={lat:f} ' +
                   '+ellps=WGS84 +datum=WGS84 +units=m +no_defs' +
                   '').format(lon=sitecoords[0], lat=sitecoords[1])

    else:
        # Set up aeqd-projection sitecoord-centered, assuming spherical earth
        # use Sphere azimuthal equidistant projection
        projstr = ('+proj=aeqd +lon_0={lon:f} lat_0={lat:f} +a={a:f} '
                   '+b={b:f} +units=m +no_defs').format(lon=sitecoords[0],
                                                        lat=sitecoords[1],
                                                        a=re, b=re)

    rad = projection.proj4_to_osr(projstr)

Returns
    -------
    r : :class:`numpy:numpy.ndarray`
        Array of xyz.shape. Contains the radial distances.
    theta: :class:`numpy:numpy.ndarray`
        Array of xyz.shape. Contains the elevation angles.
    phi : :class:`numpy:numpy.ndarray`
        Array of xyz.shape. Contains the azimuthal angles.
    """

    # get the approximate radius of the projection's ellipsoid
    # for the latitude_of_center, if no projection is given assume
    # spherical earth
    try:
        lat0 = proj.GetProjParm('latitude_of_center')
        re = projection.get_earth_radius(lat0, proj)
    except Exception:
        re = 6370040.

    # calculate xy-distance
    s = np.sqrt(np.sum(xyz[..., 0:2] ** 2, axis=-1))

    # calculate earth's arc angle
    gamma = s / (re * ke)

    # calculate elevation angle theta
    numer = np.cos(gamma) - (re * ke + alt) / (re * ke + xyz[..., 2])
    denom = np.sin(gamma)
    theta = np.arctan(numer / denom)

    # calculate radial distance r
    r = (re * ke + xyz[..., 2]) * denom / np.cos(theta)