Generate PNG/JPG maps from multiple ASCI/GeoTIFF in Python/GDAL

Question

I need to generate maps (saved in PNG, JPG or PDF) from raster files (for each file - one map) which meets following criteria; * each map has exactly the same color pallette, so that it will be easy to compare them, * I wish to add any label (for example source file name) at the top of each map, * (optionally) I wish to include legend (color pallette) * (I will be in heaven) if there will be coordinates in map included.

It is quite easy to save GeoTIFF as PNG (gdal_translate -of PNG), but I can't find a solution to set mentioned options. Is it possible to generate such maps from Python/GDAL? How? Should I call from Python QGIS map composer or there's another solution?

---

Edit

I started to write my Python script to resolve my issue. I assume that:

• I have opened a Composer View where I have all elements already added
• For each raster file, I need to:
  • remove old raster file
  • add new raster file
  • (if needed) refresh map in composer view. Legend and all other elements should probably leave the same
  • export as PDF/PNG/any other graphics format.

from PyQt4.QtGui import QImage, QPainter from PyQt4.Qt import QSize from PyQt4.QtCore import QRectF

def StringToRaster(raster):
    # This function loads raster from provided path http://gis.stackexchange.com/questions/144058/how-to-load-a-raster-layer-using-pyqgis
    if isinstance(raster,basestring):
        fileInfo = QFileInfo(raster)
        baseName = fileInfo.baseName()
        path = fileInfo.filePath()
        if (baseName and path):
            raster = QgsRasterLayer(path, baseName)
            if not raster.isValid():
                print "Layer failed to load!"
                return
        else:
            print "Unable to read basename and file path - Your string is probably invalid"
            return
    return raster

def GenerateRasterMaps(rasterList,rasterStyle):
    # I wish to use exactly the same composer view - only raster file should be changed

    myComposer=iface.activeComposers()[0]

    for rasterFile in rasterList:
        if (newRaster):
            QgsMapLayerRegistry.instance().removeMapLayer(newRaster.id())
        newRaster=StringToRaster(rasterFile)
        newRaster.loadNamedStyle(rasterStyle)
        QgsMapLayerRegistry.instance().addMapLayer(newRaster)


        # Now I should only refresh current composer view and export map as a PDF. Unfortunately, after few hours on reading (and testing): http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/composer.html and QGIS API, I still have no effect. Why it is so hard to export current Composer View to PDF/PNG. Why it is not simply one command?)

Answer 1

Since you are using QGIS, the simplest way I can think of to achieve what you need is to use Atlas Generation in Print Composer and you will have no need for any programming. Have a look at this tutorial.

Essentially, your process will be as follows:

1. Create a vector polygon file where each polygon is the extent of each raster. Do this by having all your rasters in one folder. Then use the Raster->Miscellaneous->Tile Index tool. This will create a set of polygons with an attribute called 'Location'. This attribute is the file-name and we will use this as your title later.
2. Create a Virtual Raster Catalog. This will allow you to treat all your rasters as though they were a single merged raster, while keeping the files separate. Go Raster->Miscellaneous->Build Virtual Raster (Catalog). Load and style the VRT to give you a uniform colour pallet across all rasters (alternatively, load each raster, create a style in one and then copy-paste the style to each raster by right-clicking on the layers).
3. Create a new map composer. Set up your grid, coordinates, legend, other map-furniture etc. Ensure you select the Atlas option and choose your tile index (created in step 1) as the Atlas Coverage Layer. Atlas uses the polygons as map extents and will iterate over the vector file and create a map for each polygon.
4. Add your title by adding a simple text box and use the format [%location%]. This tells Atlas to use the attribute from the vector file it is iterating over as the title of the map. In our process will be the filename of the raster (which is what I think you want), otherwise you might want to create another attribute in step 1 above, perhaps by parsing the file-paths that are set as the location attribute for each polygon by default, or if you don't have too many, you could do it by hand.
5. Finally, output the maps in one go by going Composer->Export as Image.

So now you should have a set of files (or you can export as a single multipage PDF) which are formatted maps with coordinates, a grid, north arrows and logos (if you included them), filenames as titles etc... and you should now be in heaven :)

Answer 2

You may be able to do this with mapnik (http://mapnik.org/), which has a Python API. You could create a script to loop through the TIFF files, apply styling with mapnik, and output to PNG.