def METHOD_NAME(demfile, lon, lat, outwidth=800, outheight=600):
    '''Use povray to raytrace an input Digital Elevation Model.
       demfile is a file in a format gdal can open, e.g. GeoTIFF.
       lon, lat are in degrees.
       outwidth and outheight are the size of the generated image.
    '''
    demdata = gdal.Open(demfile)
    demarray = np.array(demdata.GetRasterBand(1).ReadAsArray())
    affine_transform = affine.Affine.from_gdal(*demdata.GetGeoTransform())
    inverse_transform = ~affine_transform
    obs_x, obs_y = [ round(f) for f in inverse_transform * (lon, lat) ]
    imheight, imwidth = demarray.shape
    # In Python2, numpy can't deal with floating point indices:
    obs_ele = demarray[int(obs_x), int(obs_y)]
    lon_rad = math.radians(lon)
    lat_rat = math.radians(lat)
    # povray doesn't work if you view from ground level.
    # I don't know what the rules are or if it's better to multiply
    # or add a fudge factor.
    # multiplying by 1.15 seems to work (1.1 doesn't).
    heightfudge = 1.15
    # How big a circle, in km, can we draw around the observer
    # while staying inside the image?
    # Note that affine_transform returns degrees, not radians,
    # so that's what we need to pass to haversine_distance.
    # Edge of image north of observer:
    northlon, northlat = affine_transform * (obs_x, 0)
    northdist = haversine_distance(lon, lat, northlon, northlat)
    # Edge of image south of observer:
    southlon, southlat = affine_transform * (obs_x, imheight)
    southdist = haversine_distance(lon, lat, southlon, southlat)
    # Edge of image west of observer:
    westlon, westlat = affine_transform * (0, obs_y)
    westdist = haversine_distance(lon, lat, westlon, westlat)
    # Edge of image east of observer:
    eastlon, eastlat = affine_transform * (imwidth, obs_y)
    eastdist = haversine_distance(lon, lat, eastlon, eastlat)
    # Minimum of the distances to image edge in the four cardinal directions:
    obsradius = min(northdist, southdist, westdist, eastdist)
    print("Image size is %dx%d" % (imwidth, imheight))
    print("Observer is at (%d, %d)" % (obs_x, obs_y))
    print("Distance (km): north", northdist, "south", southdist,
          "west", westdist, "east", eastdist)
    print("Min dist:", obsradius)
    if obsradius < 3:
        print("Observer is too close to the edge")
        return
    # Loop over 8 compass points (N, NE, E etc.) and calculate the
    # pixel coordinates of the point obsradius away from the observer.
    for bearingfrac in range(8):
        bearing = bearingfrac * math.pi / 4
        bearing_deg = bearingfrac * 45
        # Find the coordinate for point with that bearing and obsradius dist:
        destlon, destlat = dest_from_bearing(lon, lat, bearing, obsradius)
        # and translate that back to pixels
        px, py = [ round(f) for f in inverse_transform * (destlon, destlat) ]
        outfilename = 'outfile%03d.png' % (bearing_deg)
        print("%3d %8.4f  %8.3f %8.3f  (%4d, %4d)" % (bearing_deg, bearing,
                                                      destlon, destlat,
                                                      px, py))
        print("%f, %f" % (float(obs_x) / imwidth, 1. - float(obs_y) / imheight))
        povfilename = '/tmp/povfile.pov'
        povfiletext = '''