From various sources I found out how to assemble a render graph to put out the depth map (scaled to the range [0, 255]):

    # Set up rendering of depth map:
    bpy.context.scene.use_nodes = True
    tree = bpy.context.scene.node_tree
    links = tree.links
    
    # clear default nodes
    for n in tree.nodes:
        tree.nodes.remove(n)
    
    # create input render layer node
    rl = tree.nodes.new('CompositorNodeRLayers')
    
    map = tree.nodes.new(type="CompositorNodeMapValue")
    # Size is chosen kind of arbitrarily, try out until you're satisfied with resulting depth map.
    map.size = [0.08]
    map.use_min = True
    map.min = [0]
    map.use_max = True
    map.max = [255]
    links.new(rl.outputs[2], map.inputs[0])
    
    invert = tree.nodes.new(type="CompositorNodeInvert")
    links.new(map.outputs[0], invert.inputs[1])
    
    depthViewer = tree.nodes.new(type="CompositorNodeViewer")
    links.new(invert.outputs[0], depthViewer.inputs[0])
    # Use alpha from input.
    links.new(rl.outputs[1], depthViewer.inputs[1])

The output can then be saved to a file with the following line

    bpy.data.images['Viewer Node'].save_render(scene.render.filepath + "_depth.png")

or create a *File Output* node:

    # create a file output node and set the path
    fileOutput = tree.nodes.new(type="CompositorNodeOutputFile")
    fileOutput.base_path = "/my_path/"
    links.new(invert.outputs[0], fileOutput.inputs[0])


You MUST start blender with the GUI for this to work. Else, the viewer will just output a black image.