Showing posts with label tip. Show all posts
Showing posts with label tip. Show all posts

Sunday, November 24, 2013

Fast tip: Drawing raster NO DATA values with MapServer

MapServer is a great way to draw web maps, using it either as a WMS server or creating image files directly.
Is the raster has a No Data value, MapServer reads it from the GDAL library and doesn't draw the pixels with the value. But what happens if you want to draw the pixels with No Data values?

Fast answer 

The solution, that I found here, is setting the LAYER property PROCESSING in this way:

PROCESSING "NODATA=OFF"
which disables the regular MapServer behaviour, which is, quite logically,  not drawing the pixels with the NODATA value.
Once this is done, just create a CLASS with the No Data Value  ¡in the EXPRESSION tag, so the pixels are coloured as you prefer.

Why did I need this stuff

Maybe this is not a common need, but when drawing weather radar images, is a good idea to indicate the areas not covered by the radar in some colour, so it's easy to distinguish between zones where there is no rain from the zones where no data is available. Like in these examples:

Catalan Weather Service radar with the range visible in gray
 The Catalan Weather Service radar in La Panadella with the range visible in gray, using the no data values
http://www.aemet.es/ca/eltiempo/observacion/radar?w=0
The Spanish Agency radar network, with the not covered areas in gray.

Strangely, finding the solution was quite difficult, so I put it in this post. Maybe somebody, someday, will find faster than me.

Some links:

Saturday, November 10, 2012

Fast tip: Changing the geotransform with GDAL python

Preparing the next post, I have found a file with a wrong geotransform, but not an easy tool to do it.
Coding it is as easy as opening the datasource with the update option and setting the new geotransform as follows:
ds = gdal.Open( fileIn, GA_Update )
ds.SetGeoTransform(geotransform)
ds = None
Where the geotransform must be a tuple like (-180.5, 1.0, 0.0, 90.5, 0.0, -1.0). Take a look to the documentation for more information.

I have created a small program to make it even easier, called changeGeotransform.py:
from osgeo import gdal
from osgeo.gdalconst import *
import sys

def changeGeotransform(fileIn, geotransform):
    """
    Takes a dataset, and changes its original geotransform to an arbitrary geotransform
    """
    ds = gdal.Open( fileIn, GA_Update )
    #Exit if there is an error opening the dataset. GDAL will output the error message
    if ds is None: 
        sys.exit()
    ds.SetGeoTransform(geotransform)
    ds = None

def helpText():
    print "To change the file Geotransform, the command is python changeGeotransform  , where the new geotransform must be like '(left_value,delta_x,rotation_x,top_value,rotation_y,delta_y)'"
    sys.exit()
def geotransformError(inputGeotransform):
    print "Error reading the geotransform: " + inputGeotransform + " it must be like (-180.5, 1.0, 0.0, 90.5, 0.0, -1.0)"
    sys.exit()


if __name__ == "__main__":
    """
    Program mainline
    """
    # Parse command line arguments.
    argv = gdal.GeneralCmdLineProcessor( sys.argv )
    if len(argv) != 3:
        helpText()
    fileIn = argv[1]
    try:
        geotransform = eval(argv[2])
    except Exception, ex:
        geotransformError(argv[2])
    if type(geotransform) != type(()):
        geotransformError(argv[2])
    changeGeotransform(fileIn, geotransform)
To execute, just type:
python changeGeotransform file_name new_geotransform
 where the new geotransform must be like:
(left_value,delta_x,rotation_x,top_value,rotation_y,delta_y)
The program basically tests if the inputs are correct and executes the code I have commented before.
Of course, only the files with drivers able to write the format can be changed. Type  gdalinfo --formats to know which drivers do you have installed.