How to use the wradlib.io function in wradlib

Returns
        -------
        ds_mem : object
            gdal.Dataset object
        """
        progress = None if (silent or isWindows) else gdal.TermProgress

        # create mem-mapped temp file dataset
        tmpfile = tempfile.NamedTemporaryFile(mode='w+b').name
        ds_out = io.gdal.gdal_create_dataset('ESRI Shapefile',
                                             os.path.join('/vsimem', tmpfile),
                                             gdal_type=gdal.OF_VECTOR)

        # create intermediate mem dataset
        ds_mem = io.gdal.gdal_create_dataset('Memory', 'out',
                                             gdal_type=gdal.OF_VECTOR)

        # get src geometry layer
        src_lyr = self.src.ds.GetLayerByName('src')
        src_lyr.ResetReading()
        src_lyr.SetSpatialFilter(None)
        geom_type = src_lyr.GetGeomType()

        # get trg geometry layer
        trg_lyr = self.trg.ds.GetLayerByName('trg')
        trg_lyr.ResetReading()
        trg_lyr.SetSpatialFilter(None)

        # buffer handling (time consuming)
        if self._buffer > 0:
            for i in range(trg_lyr.GetFeatureCount()):

silent : bool
            If True no ProgressBar is shown. Defaults to False.
        """
        silent = kwargs.pop('silent', False)
        progress = None if (silent or isWindows) else gdal.TermProgress

        layer = self.ds.GetLayer()
        layer.ResetReading()

        x_min, x_max, y_min, y_max = layer.GetExtent()

        cols = int((x_max - x_min) / pixel_size)
        rows = int((y_max - y_min) / pixel_size)

        # Todo: at the moment, always writing floats
        ds_out = io.gdal.gdal_create_dataset('MEM', '', cols, rows, 1,
                                             gdal_type=gdal.GDT_Float32)

        ds_out.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
        proj = layer.GetSpatialRef()
        if proj is None:
            proj = self._srs
        ds_out.SetProjection(proj.ExportToWkt())

        band = ds_out.GetRasterBand(1)
        band.FlushCache()
        if attr is not None:
            gdal.RasterizeLayer(ds_out, [1], layer, burn_values=[0],
                                options=["ATTRIBUTE={0}".format(attr),
                                         "ALL_TOUCHED=TRUE"],
                                callback=progress)
        else:

fcm = None
    if os.getenv('FCM_API_KEY', None):
        fcm = FCMNotification()

    # mongodb setup
    db_client = MongoClient("mongodb://mongo:27017/")
    # both will be created automatically when the first document is inserted
    db = db_client["meteocool"]
    collection = db["collection"]

    # forecast file enumeration & import
    forecast_files = sorted(glob.glob(radar_files + "/FX*_*_MF002"))
    max_ahead = 0
    forecast_maps = {}
    for f in forecast_files:
        forecast_maps[max_ahead] = wrl.io.radolan.read_radolan_composite(f)
        max_ahead = max_ahead + 5
    max_ahead -= 5
    logging.warn("Maximum forecast in minutes: %d" % max_ahead)

    # wradlib setup
    gridsize = 900

    # iterate through all db entries and push
    cursor = collection.find({})
    cnt = 0
    for document in cursor:
        try:
            doc_id = document["_id"]
            lat = document["lat"]
            lon = document["lon"]
            token = document["token"]

import wradlib
import geojson_utils
from geojson_utils import number2radius, number2degree
import sys
import math
from geojson import Point, Polygon, FeatureCollection, Feature
from PIL import Image
import colorsys

start = (46.9526, 3.5889)

img = Image.new( 'RGBA', (900,1100)) # Create a new black image
pixels = img.load() # Create the pixel map

data = wradlib.io.radolan.read_radolan_composite(sys.argv[1])
# data is mm/5m, we convert it to dbz here: https://plot.ly/~ToniBois/1783.embed
# https://www.dwd.de/DE/leistungen/radarniederschlag/rn_info/download_niederschlagsbestimmung.pdf?__blob=publicationFile&v=4

def scale(r, g, b, num, idx, a=0xFF):
    threshold = 0.2
    hls = colorsys.rgb_to_hls(r/255, g/255, b/255)
    result = colorsys.hls_to_rgb(hls[0], hls[1]+(1-hls[1]-threshold)*(1/num)*idx, hls[2])
    return (int(result[0]*255), int(result[1]*255), int(result[2]*255), a)

def dbz2color(dbz):
    factor = (dbz % 5) + 1
    if dbz >= 170:
        return (0, 0, 0, 0)
    if dbz >= 75:
        return scale(0xFE, 0xFC, 0xFD, 95, dbz % 75)
    if dbz >= 70:

"""Decodes the binary runlength coded section from DWD composite
    file and return decoded numpy array with correct shape

    Parameters
    ----------
    binarr : string
        Buffer
    attrs : dict
        Attribute dict of file header

    Returns
    -------
    arr : :func:`numpy:numpy.array`
        of decoded values
    """
    buf = io.BytesIO(binarr)

    # read and decode first line
    line = read_radolan_runlength_line(buf)
    arr = decode_radolan_runlength_line(line, attrs)

    # read following lines
    line = read_radolan_runlength_line(buf)

    while line is not None:
        dline = decode_radolan_runlength_line(line, attrs)
        arr = np.vstack((arr, dline))
        line = read_radolan_runlength_line(buf)
    # return upside down because first line read is top line
    return np.flipud(arr)

raster image with georeferencing (GeoTransform at least)
    resample : GDALResampleAlg
        If None the best algorithm is chosen based on scales.
        GRA_NearestNeighbour = 0, GRA_Bilinear = 1, GRA_Cubic = 2,
        GRA_CubicSpline = 3, GRA_Lanczos = 4, GRA_Average = 5, GRA_Mode = 6,
        GRA_Max = 8, GRA_Min = 9, GRA_Med = 10, GRA_Q1 = 11, GRA_Q3 = 12
    kwargs : keyword arguments
        passed to wradlib.io.dem.get_strm()

    Returns
    -------
    elevation : :class:`numpy:numpy.ndarray`
        Array of shape (rows, cols, 2) containing elevation
    """
    extent = get_raster_extent(dataset)
    src_ds = wradlib.io.dem.get_srtm(extent, **kwargs)

    driver = gdal.GetDriverByName('MEM')
    dst_ds = driver.CreateCopy('ds', dataset)

    if resample is None:
        src_gt = src_ds.GetGeoTransform()
        dst_gt = dst_ds.GetGeoTransform()
        src_scale = min(abs(src_gt[1]), abs(src_gt[5]))
        dst_scale = min(abs(dst_gt[1]), abs(dst_gt[5]))
        ratio = dst_scale/src_scale

        resample = gdal.GRA_Bilinear
        if ratio > 2:
            resample = gdal.GRA_Average
        if ratio < 0.5:
            resample = gdal.GRA_NearestNeighbour

# this thing is a huge mess XXX
# rename me!

import wradlib
import sys
from PIL import Image
from dbz2color import dbz2color


if __name__ == "__main__":
    start = (46.9526, 3.5889)

    img = Image.new("RGBA", (900, 1100))  # Create a new black image
    pixels = img.load()  # Create the pixel map

    data = wradlib.io.radolan.read_radolan_composite(sys.argv[1])
    # data is mm/5m, we convert it to dbz here: https://plot.ly/~ToniBois/1783.embed
    # https://www.dwd.de/DE/leistungen/radarniederschlag/rn_info/download_niederschlagsbestimmung.pdf?__blob=publicationFile&v=4

    feature_list = []
    rows = 0
    for row in reversed(data[0]):
        cols = 0
        for pixel in row:
            cols = cols + 1
            if pixel == -9999:
                # this means no data
                continue

            try:
                pixels[cols - 1, rows] = dbz2color(wradlib.trafo.rvp_to_dbz(pixel))
            except IndexError: