TekTek

I’ve just finished some work on The Weather Chaser that adds some new live satellite data for Australia & New Zealand. This data is updated hourly, usually between 15 and 25 minutes past the hour.

The first area is the addition of static ’small’ and ‘large’ sat images (both larger than the standard BoM ones) for all the different sensors (Visible, IR1, IR2, NIR & WV), as well as ‘full colour’ versions of the visible, and IR1 data that are ‘photo-realistic’. The main index page is here :

http://www.theweatherchaser.com/satellite/mtsat/

You can use the prev & next buttons to go back through the last 5 days worth of images. ‘Quick Jump’ links and an animated version, or versions of smaller images is planned.

Here is the current live ’small’ colour IR1 image (click on it to go to the large version) :

The main interactive observations map has also been updated to now have the live satellite data as a background option, including the colour versions of the Visible and IR1 sensors. This allows you to see the current satellite info overlayed with the latest radar data and observation data, and zoom/pan to the exact area of interest (and also bookmark that link).

The main page is here :

http://www.theweatherchaser.com/maps/index.html

The visible satellite layers are of the highest resolution, being similar to that of the radar data at about 1km/pixel. The other satellite layers go to about 4km/pixel.

Here is a sample direct link to some activity around Darwin with the colour visible satellite image as the background and a map and radar layer over the top. You can zoom/pan around from here. If you get to a ‘view’ you like with specific layers, zoom, position, and popup obs graph , you can bookmark the ‘Link to Current View’ link on the top right for future ease of access.

Note that the visible images stop late in the day and start early in the morning, but will show dark areas for large parts of the country from evening through to mid morning.

On a more technical note, the colour images are created using gd and python with a custom overlay method of cloud-free satellite images from 2004 with the otherwise black and white live images. With some tuning of the gamma of both images this acheives quite a nice result. The main exception is in areas of morning and evening in the visible image. In these cases the coulds tend to become a bit too transparent, and further into the darker areas where high cloud is lit by the setting sun, but the ground is not you tend to get distinctly ground-coloured clouds appearing. This is an initial version of the overlay method, further refinements may help in these areas. Nevertheless on the whole it gives a very nice result, not too far from full colour images from the Modis satellites.

Installation of djbfft didn’t ‘just work’ on my FC4 x86_64 system as per the install instructions due to an errno issue as described here for checkpassword.

At the end of the day though it is a simple change, the same as this one :

Near the top of error.h make this change :


-extern int errno;
+#include


After that it compiled perfectly, and passed all tests!

This follows on from recent posts on speed in PIL and GD in setting pixels, in particular to show a real world difference between using pure python, and a python/c solution

A common GIS task is to do simple hill shading. This is a fairly compute intensive calculation requiring specific calcultions to be done on a 3×3 pixel window for each pixel in the output image.

The best online reference I can find to the specific algorithm is this pdf. The method was published originally in Burrough, P.A., 1986. Principles of Geographical Information Systems for Land Resources Assessment. Chapter 3. Clarendon Press: Oxford.. The algorithm used here has been partially optimised relative to a completely naive solution, but could probably be optimised further.

Below is the code for this example, see the top of hillShadingPIL.py for instructions on building and running :


hillShadingPIL.py
hillShadingPIL.c
Imaging.h
ImPlatform.h


Below is a sample input image, a small Digital Elevation Model representation of Australia (0-255 ~ 0-2000m), and the result obtained with this script, both via the python method and the c method. In terms of speed, on this system here (AMDX2 4000, Fedora Core 4), the python/C based implementation overall was about 20 times faster than the pure python one, taking 0.2 version 4.3 seconds. This amounts to a very significant saving in CPU time for a relatively small investment in coding time.

Input :

Output :

See the first and second parts to this article here for a background to this post.

This test code has been further expanded. We now have reference python loop tests (ones that loop and set variables, but not images), as well as a crude ‘raw’ PIL implementation in C. What we have below are tests for different dimension images. The primary result is that the create time for PIL images seems to be longer than for gd, but that for larger images it catches up on overall process time. This is shown here by running the tests with different sizes, and can also be shown even more so by changing the code so the timing only works across the pixel setting part.

[testPixelOps] Testing with a size of 128 :
 Method          Time (s)        Times slower than fastest
+-----+         +-------+       +------------------------+
ctypes,c,GD raw 0.00014         1.0
PIL,load/c/raw  0.00041         2.9
ctypes,c,GD     0.00072         5.1
Py loop ref     0.00404         28.5
Py loop/arr ref 0.00809         56.9
PIL,load/raw    0.01829         128.7
ctypes, GD      0.08235         579.5
PIL, numpy *    0.14789         1040.7
PIL - 'point'   0.23567         1658.5

[testPixelOps] Testing with a size of 256 :
 Method          Time (s)        Times slower than fastest
+-----+         +-------+       +------------------------+
ctypes,c,GD raw 0.00038         1.0
PIL,load/c/raw  0.00099         2.6
ctypes,c,GD     0.00224         5.8
Py loop ref     0.01649         43.0
Py loop/arr ref 0.03082         80.4
PIL,load/raw    0.04195         109.5
ctypes, GD      0.17461         455.7
PIL, numpy *    0.34225         893.3
PIL - 'point'   0.84610         2208.3

[testPixelOps] Testing with a size of 512 :
 Method           Time (s)        Times slower than fastest
+-----+         +-------+       +------------------------+
ctypes,c,GD raw 0.00174         1.0
PIL,load/c/raw  0.00367         2.1
ctypes,c,GD     0.00757         4.4
Py loop ref     0.02783         16.0
Py loop/arr ref 0.05465         31.4
PIL,load/raw    0.14156         81.4
ctypes, GD      0.56828         326.8
PIL, numpy *    1.05174         604.9
PIL - 'point'   1.47204         846.6

[testPixelOps] Testing with a size of 1024 :
 Method          Time (s)        Times slower than fastest
+-----+         +-------+       +------------------------+
ctypes,c,GD raw 0.00815         1.0
PIL,load/c/raw  0.01161         1.4
ctypes,c,GD     0.04334         5.3
Py loop ref     0.11258         13.8
Py loop/arr ref 0.21994         27.0
PIL,load/raw    0.54029         66.3
ctypes, GD      2.32174         284.9
PIL, numpy *    4.31950         530.0
PIL - 'point'   5.84353         717.0

The implementation for this particular PIL test is a bit of a hack on the current PIL implementation, and will probably break for further versions of PIL, and won’t work for non ‘RGB’ or ‘RGBA’ images. It was certainly more difficult to ‘wire’ the raw-c version of the PIL test than the gd one, as it needed some looking through the PIL source, but it appears that at least for larger images it does help things along. I am surprised, in particular with exclusion of image creation costs, there is such a difference between the the two ‘top’ tests, as at the end of the day, the core loop is essentially identical C code.

Edit : Thanks to a tip from Douglas Bagnall on the python image-sig list, the reason for this discrepancy was due to the array access jumping around in memory more. In the gd example the data array is arranged with columns and rows the opposite way around. Code has been updated, and above example output re-created. Now with exclusion of image creation time things are about the same with these two cases as would be expected.

As before the source code is available in the same place, this time with an extra file, and I have included the now required PIL headers here as well for convenience. See the comments at the top of testPixelOps.py for build and running instructions :


testPixelOps.py
setPixels.c
imagecmodule.c
Imaging.h
ImPlatform.h


See the first part to this article here for a background to this post, and a third part here.

I have updated these tests to include two new methods, direct access to the pixel data in gd, and an incomplete test of the use of numpy with PIL. The conversion from a numpy array back to a PIL image doesn’t work (numpy-0.9.8, Imaging-1.1.6b2), when calling Image.fromarray I get :


...
File "/usr/local/lib/python2.4/site-packages/PIL/Image.py", line 1820, in fromarray
arr = obj.__array_interface__
AttributeError: 'numpy.ndarray' object has no attribute '__array_interface__'


However in terms of looking at this from a speed point of view it doesn’t matter, as just creation and population of the array makes it the second slowest method available. Also, the new test using raw pixel data access through gd is nearly twice as fast as using gdImageGetPixel, and hence proprtionately faster than other methods. A new indicative output is here :

> gcc -shared -lgd -fPIC -osetPixels.so setPixels.c && ./testPixelOps.py
 Method          Time (s)        Times slower than fastest
+-----+         +-------+       +------------------------+
ctypes,c,GD raw 0.00082         1.0
ctypes,c,GD     0.00177         2.2
PIL - 'load'    0.03226         39.4
ctypes, GD      0.14428         176.4
PIL, numpy *    0.26271         321.2
PIL - 'point'   0.37180         454.5

I also fixed a few bugs that affected the output, but not the timing, and added code to save the images out where possible to further check they were actually doing what was intended. So when you run the new scripts the following files should be created for comparison :


test.gd.c.gdImageSetPixel.png
test.gd.c.raw.png
test.gd.gdImageSetPixel.png
test.PIL.ImageDraw.point.png
test.PIL.Image.load.png

(They should all be 256×256 light green images)

The updated code is in the same location : testPixelOps.py and setPixels.c.

Note : There are now follow up posts to this one : here and here.

Within PIL there are two main methods for being able to set individual pixel values, via ImageDraw.point, and as of version 1.1.6, accessing the raw pixel data from Image.load as a two dimensional array. The latter being much faster than the first. It is also possible to do the same thing with GD via the ctypes library which is now standard in Python 2.5, and can be installed on earlier versions.

I have been doing a large amount of work with individual pixels, that was in turn requiring a lot of CPU time, even using the Image.load 2D Array. This had got to the point where I was running out of available CPU time on the server I had to get the work needed done. As such I thought it would be a worthwhile excercise to look at interfacing to GD via ctypes, and to compare the various speeds. One new option we have with ctypes that at least isn’t as easy to do with PIL is to create a small piece of simple C code to take over the core loops that take most of the time.

So what we have here is a script that tests four different methods of doing the same thing - setting every pixel of an image to a single colour, one at a time :

• ImageDraw.point - Known to be slow relative to the next method.
• Image.load - An improvement over the previous method.
• GD with ctypes & gdImageSetPixel.
• GD with ctypes & C based core loop.

The results of this show that using GD through ctypes and C code for the core loop can make a very significant improvement over pure python, for only a little additional complexity. Here is a sample run on an AMDx2 4000 processor running Linux Fedora Core 4 :

 Method          Time (s)        Times slower than fastest
+-----+         +-------+       +------------------------+
ctypes, c, GD   0.00181         1.0
PIL - 'load'    0.03965         21.9
ctypes, GD      0.18710         103.5
PIL - 'point'   0.45009         248.9

You can download the python file testPixelOps.py and the C code setPixels.c and try this out yourself. I would be interested to hear from anyone who can suggest a way to implement a similar loop in C for PIL to efficiently use.

At the end of the day this has resulted in saving a large number of CPU cycles, which has in turn enabled more work to be done with the same amount of hardware. It also hasn’t required a significant amount of extra work to be done from a coding point of view. I have found the approach of ‘prototyping’ in pure python, then if there is efficiency issues, moving ‘out’ the most time consuming sections where possible to C, and interfacing via ctypes to be a great way to make python code be almost as fast as a pure C application.

Note : There are now follow up posts : here and here.

Well a long overdue though very initial facelift has been finally given the main site. Primarily we now have relevant information describing services provided by Langarson including python and zope development.

A new website is now in the process of being created called The Weather Archive. The first content there is an interactive map of the seasons throughout the world. This has been created from base data from NASA’s Blue Marble project, and has been integrated using the Google Maps API. You can do most of the normal things you can do on google maps, except you get to choose the background satellite image based upon a specific month of the year. The resolution isn’t as high as the default Google Satellite data, but it does provide an interesting view into the way in which the world changes through the seasons. You can clearly see vegitation drying out, snow blanketing large areas, and swathes of green appearing. You can create and bookmark your own specific links to the page, as well as email them to yourself or friends.

>> The World Through the Seasons

I have been doing a lot of work in Javascript in particular with Firefox of late (DHTML & extensions), here are a few things I have learnt along the way in terms of being able to effectively debug and develop this code. In particular I am new to the extension development world which has it’s own set of differences from plain javascript in a web page.

Where is that javascript error?

After much frustration I stumbled across the reason why the javascript console wasn’t showing all my javascript errors. To get errors within all the code, or within your extensions showing inside Firefox you need set ‘javascript.options.showInConsole’ to true. Put in ‘about:config’ in the location bar, find that entry, toggle it to true, and restart your browser. Now javascript errors that were being swallowed can be seen on the javascript: console. The other little thing I only just noticed was a context menu (right click) on the console allows you to put the most recent at the top, I had always been frustrated having to scroll to the bottom if there were errors reported in other parts of the browser I had to scroll over.

The new javascript console
I have also just started using a new javascript console - Console². This is quite promising, as one big gripe about the normal console is the fact you can’t filter errors by their type. Anyway, have a look at this thread on mozillazine for some more information, and updates.

Restart, restart, restart

Perhaps the other big help is the ability to run multiple instances of the browser so you can easily restart it (or crash it!) when developing for extensions without closing a browser window you may have a number of reference/research open in. Create a new profile say ‘test’ then on windows at least, create a new ‘.bat’ file with content along the lines of :

@echo off
set MOZ_NO_REMOTE=1
start “” “C:Program Files\Firefox 1.5\firefox.exe” -p “test”
set MOZ_NO_REMOTE=0

Put that somewhere handy on your desktop/quick start bar/start menu, and you can easily develop with two browsers. Add more profiles if you wanted more browsers. The other good thing about this is you can have this profile configured with minimal extensions to perhaps help with startup time, and even make the default page ‘about:blank’ as well.

A suggestion is to change the theme in this test profile to be visually very different so that you know which one it is, and don’t start closing the wrong one. Another suggestion is to setup a different icon as well so you can easily see from the taskbar which is which. To do this create this folder if it doesn’t exist :

C:\Program Files\Mozilla Firefox 1.5\chrome\icons\default

And put your new icon in there called main-window.ico. This does mean you need two installs of firefox though if you want different icons for the two profiles which is a bit of a pain. It would be nice if this could be overridden in the user settings somehow, it perhaps can, but I haven’t investigated that as yet.

Venkman the Javascript debugger

For normal javascript debugging this is a great tool. Here is the site for 1.5+ : http://getahead.ltd.uk/ajax/venkman

In order to get it to debug extensions though you need to do deselect the default selected ‘Debug->Exclude Browser Files’. Similar to the javascript errors not showing up I was wondering how to get the ‘browser’ files in here, and it turned out to be quite obvious in the end. One thing that will become more useful with this turned on though is the search facility for the loaded scripts at the top of the list as all of a sudden you will see a lot more files in this list.

Other resources
There are a few other useful notes on the main mozilla developer site that are worth looking at :
http://developer.mozilla.org/en/docs/Building_an_Extension#Debugging_Extensions

I have been working on a couple of google maps projects, and in particular ones that overlay a custom layer of live weather radar data. Currently this data is being sourced from the Australian BoM, and the NOAA in America. As such coverage extends over Australia, and the US. Links to each of these projects :

• Interactive Australian Observations Map
• Interactive US Observations Map

Creating this custom overlay was greatly assisted by the resources such as this one at mapki. At the end of the day however there were a number of challenges left to get this thing to fly. These included :

• Obtaining all the required data to create the overlay :
• Both Australia and the US have a lot of radar stations (total around 200).
• Data for each station can consist of a number of files.
• Maintaining the current set of these consumes what at least currently would be considered reasonable amount of bandwidth (Particularly in Australia).
• Efficiently creating live overlay tiles from data that changes every 5-10 minutes :
• Converting the original data into a single data set with the right projection (Mercator) required quite a bit of work.
• Then turning this base data into the tiles required more custom development, and a reasonable amount of server resources. All existing tools I could find were either not suitable for running on a *nix based system, or they were not going to be as efficient as something written in plain C.
• The first effort was completely written in python using PIL on the server side, but it quickly became apparent that it was not going to be quick enough to create the final product in a timely enough fashion, whilst leaving some spare CPU cycles on the server. To resolve this problem core image processing routines were re-written in straight C and used the gd image library and connected to the python wrapper scripts using SWIG. The result was a 5-20 times speedup of the code, and meant the server had plenty of ’spare’ CPU cycles.
• Creating the tiles in just the right format :
• The main little niggle here was the fact that images produced out of PIL at least (not tested with gd) where the 0 palette entry was black, and set as the transparent value would come up as white in Internet Explorer (IE). I didn’t fiddle too much with this, as I found if I used another colour not used by the radar intensities that it would work fine in IE, and still worked fine in Firefox. So I don’t know now wether it was a problem with PIL, IE, or both.
• One other related thing here was the importance of maintaining properly paletted images through out the processing chain. In general using ‘True Colour’ (24 bit) images used a lot more memory, and also caused problems when converting back to palettes. This included undesired dithering of GIF images when saved via PIL without first being specially converted to a palleted image. ie :

image = image.convert('P',palette=Image.ADAPTIVE,colors=256, dither=Image.NONE)

image.save(fp, 'GIF')

Many of the points above could have a whole article, so what I hope to do is to elaborate on some of these and publish some code/utilities related to them over the coming weeks, particularly as I do further refinement of them