Raster calculations with GDAL and numpy: calculating GFS wid speed

Performing raster calculations is a frequent need. GDAL python bindings come with a utility program,  gdal_calc.py, that does that, but many times is useful to code them.

In this example, we will be able to calculate the wind speed field from the GFS model, from the given u and v speed.

 The calculated wind field. America is easily seen at the right side of the image

The data

In this example, we will use the GFS global model data, which is available at their web site. Just navigate to a date that pleases you and download a .grb file.

Files are in grib format, with a little problem when trying to view them in QGis: Coordinates are from 0 to 360º in longitude, while the usual thing is to put them from -180º to 180º. Everything will work without changing this, but you can do it after taking a look to this post.

Inside the grib file, there are many bands (hundreds). You can see their definition using the gdal command gdalifo, or the file with the extension *inv at the same place you downloaded the grib file. There are many wind fields, that correspond to different altitudes, all them indicated as u (east) and v components (north). Just take note of the band number of both at the desired level.

NumPy

The best of using the python gdal bindings, compared from those of the other programs, is that converting a raster file to a NumPy matrix is now very easy. 

NumPy is a python extension that allows to perform matrix calculations in a really easy and efficient way. It should be installed in your computer if GDAL python is installed.

Once you have two matrices A and B in NumPy, adding element by element to a C matrix would be as easy as:

C = A + B

and multiplying:

C = A * B

Note that is not a matrix multiplication (more info here), but this is what we usually need when working with raster images.

The code

The wind speed is the modulus of the two components, so 

spd = sqrt(u^2+v^2)

The "traditional way" to calculate it would be opening two bands and loop all the elements to create the new band. Instead, this can be done:

from osgeo import gdal
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *

fileName = "gfs.tiff"
bandNum1 = 199
bandNum2 = 200

outFile = "out.tiff"

#Open the dataset
ds1 = gdal.Open(fileName, GA_ReadOnly )
band1 = ds1.GetRasterBand(bandNum1)
band2 = ds1.GetRasterBand(bandNum2)

#Read the data into numpy arrays
data1 = BandReadAsArray(band1)
data2 = BandReadAsArray(band2)

#The actual calculation
dataOut = numpy.sqrt(data1*data1+data2*data2)

#Write the out file
driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create("out.tiff", ds1.RasterXSize, ds1.RasterYSize, 1, band1.DataType)
CopyDatasetInfo(ds1,dsOut)
bandOut=dsOut.GetRasterBand(1)
BandWriteArray(bandOut, dataOut)

#Close the datasets
band1 = None
band2 = None
ds1 = None
bandOut = None
dsOut = None

Note that:

1. gdalnumeric is used, because it's the easiest way to dump a raster into a NumPy matrix using BandReadAsArray and the inverse function with BandWriteArray.
2. In this case, we are reading the bands 199 and 200, but you should check which bands are the ones to read in the file you download.
3. CopyDatasetInfo will copy all the input file metadata into the output file. This includes the geotransform, the projection, the GCPs (there aren't in this case) and the other metadata.
4. Don't forget to set the bands and datasets to None at the end to free them. In longer scripts this may be critical.

gdal_calc.py

If you want to use the command line function to do the same, you can try with:

gdal_calc.py -A gfs.tiff --A_band 199 -B gfs.tiff --B_band 200 --outfile out.tiff --calc="sqrt(A*A+B*B)"

It gives the same result.

kartograph tutorial II: going interactive

In my last post I showed how to create a cool SVG map using Kartograph, to be opened and edited with an Inkscape-like program.
The elections were celebrated and CiU won again, but with less votes. I could either update the map or play with the same data. The fun with Kartograph is using the generated SVG to create interactive web pages using Kartograph.js, so this is what I'll do.
The web page will be able to:

• Generate the map not only for the first position, but also for the others
• Query every municipality to know which were the results there.

As usual, the code is available to download, and at the working web page too.

This second example has a continuation post: Kartograph tutorial III: Symbols

 

Changes to the SVG

In the other post, we used the SVG only to draw, but now we need to add data to it, so we can query it.

So change the municipalities layer at the file elections.json from:

   "municipalities":{
       "src": "./mun_out.shp"
   }

to:

   "municipalities":{
       "src": "./mun_out.shp",
       "attributes": {"winner": "Winner","ciu": "CiU","psc":"PSC-PSOE","erc":"ERC","pp":"PP","icv":"ICV-EUiA","cs":"C's","name":"Name"}
   }

This will add the attributes set in the shapefile geometries to the SVG geometries.
To generate the SVG again, just type:
kartograph -o elections.svg elections.json

Viewing the map in a web page

Even though is quite simple, I haven't found any really basic working example, so here it is:

<html>
    <head>
        <script src="jquery-1.8.3.min.js"></script>
        <script src="raphael-min.js"></script>
        <script src="kartograph.min.js"></script>
 <script language="JavaScript">
function loadMap(){
var map = Kartograph.map('#map', 600, 0);

map.loadMap('elections.svg', function() {

    map.addLayer('background');
    map.addLayer('world');
    map.addLayer('depth');
    map.addLayer('trees');
    map.addLayer('crops');
    map.addLayer('municipalities');
     
}, { padding: -30 })
}
</script>
    </head>
    <body onLoad="loadMap()">

        <div id="map"></div>
    </body>
</html>

1. Include jquery, raphael and kartograph
2. The map must be created when the page loads, so the function loadMap() is called at the body tag , using the onLoad event.
3. Kartograph.map('#map', 600, 0) creates the map at the div tag with the id map. X size will be 600 px, and the aspect ratio will be respected, because the y size is 0. When an other y is set, strange things happen with the png borders.
4. loadMap draws the map from the file elections.svg. Then, executes a function. In this case, the function loads all the layers one by one, using the method addLayer(LayerNameInTheSVG)
5. padding: -30 removes the strange lines from the vectors in the borders by padding the image. Try without using it to understand what happens.

This renders the SVG without any styling or interaction. But is the basic way to start any Kartograph.js project.

CSS Styling

There are two ways to add styles to the map. We will use both, to show them.

The first one, is to use a separate css file, as in the last post.

Adding the css file is quite tricky, because of the IExplorer compatibility. Now, the loadMap call would look like:

map.loadMap('elections.svg', function() {
 map.loadCSS('elections.css', function() {
  map.addLayer('world');

   ...
 });
})

So, the only thing that the loadMap method does, is calling the loadCss method, which will do all the stuff.

The css file is also different than the one used before. To style the layer world, the css section will be:

#map svg path.world {
stroke: #fff;
fill: #cccccf;
fill-opacity: .7;
}

Tip1: Be careful, Chrome doesn't reload your css even though you use ctrl-F5!
Tip2: Do not separate the CSS elements with comas (look at the example). With the python example, this didn't matter.
The other way is using the style method from the layer, as we will see in the final part of the post.

Adding tooltips

If we want the web page to show a tooltip with the result when the mouse is over a municipality, just add a tooltip to the desired layer (municipalities, in this case). The function in loadCSS will be now:

map.loadCSS('elections.css', function() {
 ...
 map.addLayer('municipalities');

$.fn.qtip.defaults.style.classes = 'qtip-light';

 map.getLayer('municipalities').tooltips(function(data) {

return [data.name, 'CiU: <b>' + data.ciu + '</b><br/>PSC-PSOE: <b>'+ data.psc + '</b><br/>ERC: <b>' + data.erc+'</b><br/>PP: <b>'+data.pp+'</b><br/>ICV-EUA: <b>'+data.icv+'</b><br/>C\'s: <b>'+data.cs+'</b>'];
});
});

A function is called as a parameter to the tooltips method of the layer municipalities. This function returns what we want to be shown in the tooltip. The data parameter has all the attributes in the SVG geometries. In this case, a list with the results is shown.
The other line in the code just sets the style of the tooltip. The official documentation has more information about the options.
Also, the tooltips styles and javascript must be included in the header, adding:

<script src="jquery.qtip.js"></script>
 <link rel="stylesheet" type="text/css" href="jquery.qtip.css"/>

Setting the colors in the map

Finally, the map will be colored. By default, depending on the winner party, but with a combo box to choose the position.

So first, a function is created to sort the results in every municipality, and set the style with a different color for every party:

function setColors(value){
 map.getLayer('municipalities').style('fill', function(data) {

    var positions = new Array();
    positions[0] = parseInt(data.ciu);
    positions[1] = parseInt(data.psc);
    positions[2] = parseInt(data.erc);
    positions[3] = parseInt(data.pp);
    positions[4] = parseInt(data.icv);
    positions[5] = parseInt(data.cs);

    positions.sort(function(a,b) { return b-a; });

    var bgColor;
    if (positions[value-1] == parseInt(data.ciu)){
        bgColor = "#99edff";
    } else if (positions[value-1] == parseInt(data.psc)){
        bgColor = "#ff9999";
    } else if (positions[value-1] == parseInt(data.erc)){
        bgColor = "#EDE61A";
    } else if (positions[value-1] == parseInt(data.pp)){
        bgColor = "#005aff";
    } else if (positions[value-1] == parseInt(data.icv)){
        bgColor = "#74ff74";
    } else if (positions[value-1] == parseInt(data.cs)){
        bgColor = "#ff7b00";
    }
    
    
    return bgColor;
 });
}

First, we set the style to the layer with the method style. the first argument is the property to change, fill in this case. The second parameter is the value to set. In this case, it will depend of the function that chooses the position. Things to look at:

• value is the position to show (from 1 to 5)
• To sort the positions in the array in descending order, use the sort() method, with a function to order the elements as numbers. See why here.
• Then, check which party corresponds to the number of votes at the desired position, and set the color.

This function will be called at the beginning to show the first position by default:

setColors(1);

and a select element is created to select the position and change the colors when the option is changed:

<select onChange="changePos(this)">

<option value="1">1st position</option>

<option value="2">2nd position</option>
<option value="3">3rd position</option>
<option value="4">4th position</option>
<option value="5">5th position</option>
</select>

An auxiliar function is created to get the option value:

function changePos(sel){
 var value = sel.options[sel.selectedIndex].value;
 setColors(value);
}

That's all!
The best is to download the code or look at the working web.
Of course, some styling is still needed to do a deliverable web page, but the main idea is already done.

Kartograph tutorial: Electoral map

Tomorrow I will be all the day at the polling station, since I've been chosen as a member.
The last two weeks I was playing with the amazing Kartograph software, so it was a good moment to experiment with electoral maps (the first time for me).

 In this example, I will explain step by step how to create the map above. It's quite similar to this tutorial, but I want to continue in a new post, going interactive.
Kartograph creates vector SVG images, which you can edit later with Inkscape or Illustrator, so gives much more flexibility than systems generating PNG like files, much more difficult to modify.
As in all the posts, you can get all the files used in the example.

This example has two continuation posts: 

• Kartograph tutorial II: going interactive
• Kartograph tutorial III: Symbols

Getting the data

As usual, getting the data is not so easy. The Catalan government has a good open data web site, where I found:

• A file with a really descriptive name, bm50mv33sh1fc1r170.zip, with all the administrative boundaries (provinces, comarques, and municipalities).
• Lots of files with election results. I choose the 2010 elections, since they where to the Catalan parliament, like the ones tomorrow. As you can see on the map, the party CiU won with a very big majority, so the map is not as interesting as it could be.

I have used the municipalities to draw the map because the result is more diverse than using bigger zones. Actually, the real  constituency is the province, but CiU won everywhere, and a plain blue map is quite boring.
So I've had to join the two files to get one file with the geometries and the results. The process is quite long and dirty (why didn't they use an id? I had to join with the names), so I won't explain how to do it, but put the result at the data files. You can find this file here.

Then, to decorate the map, I used the following files

• World boundaries from Natural Earth (ne_10m_admin_0_countries.zip), to draw the coast line outside Catalonia
• From VMAP, the layers Trees, Crops, and DepthContours, to decorate the map outside the electoral  constituencies.

Since the layers are worldwide, so very big, I have used these ogr commands to clip:
 ogr2ogr -clipsrc -3 37 4 44 Trees2.shp Trees.shp
and to simplify:
 ogr2ogr -simplify 10 munis.shp bm50mv33sh1fpm1r170.shp
Doing so, the time to generate the map is divided by five or more.

Installing Kartograph

Since we only need kartograph.py for this tutorial,  first, download it from the github page clicking at the zip icon.

In a linux system, uncompress and execute 

python setup install

as a super user.

That's all, if you have the GDAL python bindings installed.

Creating the map

To create a map with Kartograph, you will need a configuration file in JSON format, which will have three basic sections:

Projection

To set the projection, there used to be a web page named Visual map configurator, that doesn't work any more. But don't worry, you can use the Map Projections page. Just choose the projection that fits you more, change the parameters and click the gear icon:

A dialog will open, and the lines that are interesting in this case are, in the image example, like:

         "proj": "sinusoidal",
        "lon0": 20

This will be translated in our json file as:
     "proj": {
            "id": "sinusoidal",
           "lon0": 20
     }

 Bounds:

The part of the world we want to represent is set here. It's quite well explained at the documentation, but it can be a bit confusing, and not all the options work with all the projections.

In our example, I have used:

  "bounds": {
    "mode": "bbox",
    "data": [-0, 40, 4, 43],
    "crop": [-3, 37, 5, 44]
  } 

• mode: How the bounds are expressed. BBOX is the basic option, but you can also set it defining the points you want to enter in the map, or even the features in a layer. If the layers are in different projections, other modes can be a little tricky.
• data: In our case, the bounding box. In other modes, the layer name, the points, or whatever.
• crop: Is an optional label. Our svg will be clipped at the bounds set at data, but all the data in the files will be processed. If the files include all the world, this takes a long time, and generates much bigger SVG outputs. With crop, only the features inside the BBOX will be included.

Layers:

As the name suggests, the layers to include. 

The shapefiles are added as:

   "municipalities":{
       "src": "./mun_out.shp"
   }

There are also two special layer, graticule and sea. The first draws the meridians ans parallels, while the second does nothing more than giving a feature to draw the background:

   "background": {"special": "sea"},
   "graticule":{ "special": "graticule", "latitudes": 1, "longitudes": 1}

All the layers  will be drawn in the order indicated at the json file, so this must be well chosen to select which layer hides what.

Styling

This is the nice part. Without styling, the SVG can be used directly with Inkscape or Kartograph.js, but is possible to generate styled maps directly with kartograph.py.
You can give the style either in the json file or in a separate css file, which seems cleaner. The names given to the layer are the ones to be used in the css as the id. So to give a style to the municipalities layer, add

#municipalities {
 fill: #FFF;
 stroke: #882222;
 stroke-width: 0.5px;
 stroke-opacity: 0.4;
}

The general options are at the documentation again. CSS for SVG is a little different from the one used in traditional html.
Since we want to paint the municipalities in a different color depending of the party who won the elections, we will use filters, like this one:

#municipalities[Winner=CiU]{
 fill: #99edff;
}

It would be nice to compare different fields i.e. CiU > PSOE, but this is not possible (at least, I haven't found how to do it), so I had to calculate the winner and put it in a field (called Winner, as you can see in the example)

Drawing

There are two options to draw the map. A command line program is installed with the setup, called kartograph. 

To draw the styled map, just type

   kartograph elections.json --style elections.css -o elections.svg

But you can also include all this in a python program, so could generate the data and then the map. In our case, the code would be

from kartograph import Kartograph
from kartograph.options import read_map_descriptor
import sys
K = Kartograph()
css = open("elections.css").read()
cfg = read_map_descriptor(open("elections.json"))
K.generate(cfg, outfile='elections.svg', format='svg', stylesheet=css) 

Finally, I edited the svg file with Inkscape to put the titles and legend. Is just to show that the idea is generating a base svg and from there, draw the pretty final map.

Configuration files

To draw the map in the example, I have used the following files:

elections.json

{
"proj": {
        "id": "sinusoidal",
        "lon0": 20
  },
   "layers": {
   "background": {"special": "sea"},
   "graticule":{ "special": "graticule", "latitudes": 1, "longitudes": 1, "styles": { "stroke-width": "0.3px" } },
    "world":{
       "src": "data/ne_10m_admin_0_countries2.shp"
   }, 
   "trees":{
      "src": "data/Trees2.shp",
      "simplify": true
   },
   "crops":{
      "src": "data/Crops2.shp",
      "simplify": true
   },
   "depth": {
       "src": "data/DepthContours2.shp",
       "simplify": true
   },
   "municipalities":{
       "src": "./mun_out.shp"
   }
   },
  "bounds": {
    "mode": "bbox",
    "data": [-0, 40, 4, 43],
    "crop": [-3, 37, 5, 44]
  }
}

elections.css

#background {
 fill: #e8f9fb;
 stroke: none;
},
#world {
 fill: #f5f3f2;
 stroke: none;
},
#graticule {
 stroke-width: 0.3px;
},
#municipalities {
 fill: #FFF;
 stroke: #882222;
 stroke-width: 0.5px;
 stroke-opacity: 0.4;
},
#municipalities-label {
 font-family: Arial;
 font-size: 13px;
},
#municipalities[Winner=CiU]{
 fill: #99edff;
},
#municipalities[Winner=PSC-PSOE]{
 fill: #ff9999;
},
#municipalities[Winner=ERC]{
 fill: #EDE61A;
},
#depth {
 stroke: #223366;
 stroke-width: 0.5px;
 stroke-opacity: 0.4;
},
#trees {
  fill: #d2f8c0;
  stroke: none;
},
#crops {
  fill: #fcf8d8;
  stroke: none;
}

What's next

If I have time, I'll try my first Kartograph.js example. From the svg generated, is possible to create cool interactive maps.

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.

IRIS driver for GDAL

IRIS is a software by VAISALA (formerly by Sigmet) to manage and visualize meteorological radar data. It has its own data format, which is pretty well documented at their web page.

Since I've worked for a long time at the SMC, I've had to work with these kind of stuff many times, developing scripts and programs to change its format so it was possible to draw the radar data.

Besides, I started using GDAL quite a lot, and wondered why it wasn't possible to open IRIS data files directly. Since nobody did a driver for it, I decide to do it myself after reading the GDAL driver tutorial, which seemed easier than what it actually was to me.

But finally I've managed to create the driver. When compiled into the GDAL distribution, the driver allows to change IRIS data format to GTiff or any other using gdal_translate, opening the file using QGIS, serving the data with Mapserver, etc. So your life will become much easier if you have to work with IRIS data files, without having to code the format translation yourself.

This first version comes with a CAPPI example file to test it,  and supports the following data types:

• PPI (reflectivity and speed): Plan position indicator
• CAPPI: Constant Altitude Plan position indicator
• RAIN1: Hourly rainfall accumulation
• RAINN: N-Hour rainfall accumulation
• TOPS: Height for selectable dBZ contour
• VIL: Vertically integrated liquid for selected layer
• MAX: Column Max Z WF W/NS Sections

I've submitted the driver at the GDAL trac so they can include it to the software if they find it convenient: https://trac.osgeo.org/gdal/ticket/4854 and commented it at their mailing list.
 But if you are impatient (or the driver isn't finally included), it is possible to add it to your GDAL after re-compiling it. The instructions:

• Download the latest GDAL source and the driver source code.
• Copy the iris folder into the frmts folder in your GDAL distribution.
• Then following the steps at the GDAL driver tutorial:
•  Edit the file gcore/gdal_frmts.h  and add the driver at the end of the list, which will look like:
  void CPL_DLL GDALRegister_IRIS(void);
• Edit the file frmts/gdalallregister.cpp and register the driver adding a piece of code:
  #ifdef FRMT_iris
      GDALRegister_IRIS();
  #endif
• Edit the files GDALmake.opt and GDALmake.opt.in and look for the line containing GDAL_FORMATS =, adding the word iris to the formats list in the variable.
• Edit frmts/makefile.vc and change the EXTRAFLAGS = list adding
  -DFRMT_iris
•  Once all this is done, just compile GDAL as usual:
  ./configure (add the flags you need here, IRIS format does not need any extra flag)
  make
  make install

 And that's all. As you can see, it would be much easier if the driver was included in the GDAL distribution!
If you have any comment or problem, please write it in the comments.

The latest SMC radar image processed from the IRIS data format