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

Keyword Arguments
    -----------------
    source_srs : osr.SpatialReference
        Source Projection

    Returns
    -------
    geom : ogr.Geometry
        Transformed Geometry
    """
    gsrs = geom.GetSpatialReference()
    srs = kwargs.get('source_srs', gsrs)

    # srs is None assume wgs84 lonlat, but warn too
    if srs is None:
        srs = projection.get_default_projection()
        warnings.warn("geometry without spatial reference - assuming wgs84")

    # transform if not the same spatial reference system
    if not srs.IsSame(dest_srs):
        if gsrs is None:
            geom.AssignSpatialReference(srs)
            gsrs = geom.GetSpatialReference()
        if gdal.VersionInfo()[0] >= '3':
            dest_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
            gsrs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
        geom.TransformTo(dest_srs)

    return geom

phi_m = np.radians(lat)
        lam_0 = 10
        lam_m = lon
        lam = np.radians(lam_m - lam_0)
        er = 6370.040
        m_phi = (1 + np.sin(phi_0)) / (1 + np.sin(phi_m))
        x = er * m_phi * np.cos(phi_m) * np.sin(lam)
        y = - er * m_phi * np.cos(phi_m) * np.cos(lam)
    else:
        # create radolan projection osr object
        proj_stereo = projection.create_osr("dwd-radolan")

        # create wgs84 projection osr object
        proj_wgs = projection.get_default_projection()

        x, y = projection.reproject(lon, lat, projection_source=proj_wgs,
                                    projection_target=proj_stereo)

    return x, y

if trig:
        # calculation of x_0 and y_0 coordinates of radolan grid
        # as described in the format description
        phi_0 = np.radians(60)
        phi_m = np.radians(lat)
        lam_0 = 10
        lam_m = lon
        lam = np.radians(lam_m - lam_0)
        er = 6370.040
        m_phi = (1 + np.sin(phi_0)) / (1 + np.sin(phi_m))
        x = er * m_phi * np.cos(phi_m) * np.sin(lam)
        y = - er * m_phi * np.cos(phi_m) * np.cos(lam)
    else:
        # create radolan projection osr object
        proj_stereo = projection.create_osr("dwd-radolan")

        # create wgs84 projection osr object
        proj_wgs = projection.get_default_projection()

        x, y = projection.reproject(lon, lat, projection_source=proj_wgs,
                                    projection_target=proj_stereo)

    return x, y

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)

    r = np.asanyarray(r)
    theta = np.asanyarray(theta)
    phi = np.asanyarray(phi)

    if r.ndim:
        r = r.reshape((1,) * (3 - r.ndim) + r.shape)

    if phi.ndim:
        phi = phi.reshape((1,) + phi.shape + (1,) * (2 - phi.ndim))

    if not theta.ndim:
        theta = np.broadcast_to(theta, phi.shape)

    dims = 3
    if not strict_dims: