Applying color style to GeoTIFF layer using Mapnik and Python

I am trying to apply a color style to a GeoTIFF using Mapnik.

The code I’m using is:

import mapnik

map = mapnik.Map(800,600)
map.srs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

style = mapnik.Style()
rule = mapnik.Rule()
sym = mapnik.RasterSymbolizer()
 
rc = mapnik.RasterColorizer()
rc.colorizer = mapnik.RasterColorizer(mapnik.COLORIZER_LINEAR, mapnik.Color("white"))
rc.colorizer.add_stop(-20, mapnik.Color("gray"))
rc.colorizer.add_stop(0, mapnik.Color("yellow"))
rc.colorizer.add_stop(10, mapnik.Color("red"))
rc.colorizer.add_stop(15, mapnik.Color("green"))
rc.colorizer.add_stop(25, mapnik.Color("blue"))


sym.colorizer = rc
rule.symbols.append(sym)
style.rules.append(rule)

map.append_style('style', style)

layer = mapnik.Layer('My data', map.srs)
layer.datasource = mapnik.Gdal(file='input.tif')

layer.styles.append('style')

map.layers.append(layer)


mapnik.render_to_file(map, 'map-output.png', 'png')

The GeoTIFF gdalinfo output is:

$gdalinfo input.tif
Driver: GTiff/GeoTIFF
Files: input.tif
Size is 1414, 631
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-125.229333333000000,50.146000000000001)
Pixel Size = (0.041700000000000,-0.041700000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-125.2293333,  50.1460000) (125d13'45.60"W, 50d 8'45.60"N)
Lower Left  (-125.2293333,  23.8333000) (125d13'45.60"W, 23d49'59.88"N)
Upper Right ( -66.2655333,  50.1460000) ( 66d15'55.92"W, 50d 8'45.60"N)
Lower Right ( -66.2655333,  23.8333000) ( 66d15'55.92"W, 23d49'59.88"N)
Center      ( -95.7474333,  36.9896500) ( 95d44'50.76"W, 36d59'22.74"N)
Band 1 Block=1414x1 Type=Float32, ColorInterp=Gray
  Min=-23.065 Max=24.441
  Minimum=-23.065, Maximum=24.441, Mean=0.561, StdDev=6.495
  NoData Value=-3.39999999999999996e+38
  Metadata:
    STATISTICS_MAXIMUM=24.44076094814
    STATISTICS_MEAN=0.5611245629664
    STATISTICS_MINIMUM=-23.065323150453
    STATISTICS_STDDEV=6.4950334131849

The problem is that the PNG result does not show the colors I’m applying. It shows only a empty PNG image.

layer.datasource = mapnik.Gdal(file='input.tif', band=1)