Read online: https://github.com/codingisacopingstrategy/raduga-server/blob/master/RAINBOW-ALGORITHM.md
Rainbows need rain and sun. So we need to figure out where is the sun, and where is the rain. Firstly we figure out where it rains.
Even if there are many ways in which meteorologists measure and predict rain, including rain gauges and weather balloons, Radars are one of the most popular means of following the development of rainclouds. In Holland, this fact has entered public consciousness through a popular website called Buienradar. Specific areas of the earth are well covered by radar. Free sources for such data are available in the US (example), Scandinavia & the Baltic states, the UK. Such a radar network doesn’t exist for Russia: the available Radar sources cover the main airports. The alternative is to use Satellite data. The first source we tried is EUMETSAT, the European Meteorological Satellites Agency. Their Multi-Sensor Precipitation Estimate (available as a Mercator projection as well) proposes ‘Operational weather forecasting in areas with poor or no radar coverage’. The EUMETSAT image doesn’t cover Russia really well though. The European satellites are in a geostationary orbit that always gives them a view of the same slice of earth, and funnily enough that slice does not cover Russia. We wondered if there are historic, old NATO style political reasons for that?
A telephone adventure calling everyone from the KNMI to Meteo-consult to the Finnish Meteorological Agency to the World Meteorological Organization (Geneva), and they didn’t know where to get data on the Russian clouds. With the help of a translator we contacted the Russian Weather agency. This helped us forward, as they pointed us to their data: Global medium-range forecast fields in GRIB format. However, until know we had imagined to work with images This meant using GRIB data—a data format that while a standard in the world of Meteorological analyses, we did not know how to use directly. Luckily we found a great example in Cameron Beccario’s visualisation called earth. For this animated map of global wind, weather, and ocean conditions, Beccario not only released the source code, but also added a detailed README explaining how he gets the source GRIB data and how to read it with a utility developed by him (GRIB2JSON).
We created a script to download the necessary GRIB files. The details can be found in fetch.py. There are some 300 different tables in the GFS data, so we ask it to filter for Precipitable Water (PWAT), the information that interests us. The GFS data is produced every six hours. It contains information about the current weather situation and for 3, 6, 9, 12 etc. hours into the future. Because the GFS is not immediately available (in general several hours after the time indicated as ‘now’), we use the predictions of 6 hours and 9 hours into the future. This way we have a time point for every 3 hours.
Now that we have the GRIB data about precipitation, we can get to the next step: processing this data. This happens in the file water.py. First we convert the GRIB to a JSON file, that we can easily read in with python. As part of this JSON, we have a long list of values that represent the amount of precipitable water for each coordinate on a 0.5 degree grid spanning the earth. This data, we convert into a form that is more easily manipulable still: an image.
Since the GRIB data already provides us with information about the coordinate system, so we can use this to map to an image in a straightforward fashion. We only have to transform the values so they map to the 0-255 range of a grayscale image. From here on we treat the data as an image. An image is just a two-dimensional array, but thinking about the data as an image as opposed to a series of numbers makes it easy to conceptualise the transformations and to show them.
To see a rainbow, one needs to be standing with ones back to the sun, looking in the direction of the raincloud. The sun can not be too high. Water always refracts light in the same way, the angle between the incoming sunlight and the rainbow being 180 minus 42 degrees. Once the sun gets to high (above 42 degrees) the rainbow disappears behind the horizon. That is why you see more rainbows in spring and autumn, as the sun stays lower throughout the day.
For our given point in time, we need to determine where the sun is between 0 en 42 degrees. The first step is to calculate the solar altitudes. This is different for each date, and it is not the same each year. We use the Pysolar module to do this. We re-use the coordinate system of the precipitation data, we loop over every latitude and longitude combination and ask Pysolar to GetAltitudeFast(latitude, longitude, DATE).
The approach might not be the fastest, but it allows us to easily find the location of the sun (there where the sun is straight overhead the earth): it is the pixel with the highest solar altitude.
We also create a mask, that blacks out all parts where the solar altitude is not between 0 and 42 degrees.
Rainbows can appear at the edge of rainclouds. But they can only appear if the rainclouds find themselves in the opposite direction of the sun.
To find these edges, we’ve used the following approach:
- Take the rain-information generated in step 1. Contrast and treshold it as to generate a mask:
- Copy the mask, and slightly extrude it with the position of the sun (found in step 2) as the center point
- Subtract the original clouds from the extrusion
What we have left is the possible position of the rainbows. Then
- Subtract the mask of sun-positions generated in step 2
This gives us all possible rainbow-areas on earth. For the purposes of Радуга, we:
- Mask out Russia
Because at this point we are manipulating images, we can use the Python Imaging Library for this purpose. Except for the extrusion: I couldn’t find out how to use PIL for this so we do it with the classic ImageMagick command line software.
For the areas where the sun is shining, we find the edges of the rainclouds from the center of the sun.
To make the rainbows visible, we have to encode the pixels that represent rainbow areas so that they can be displayed on an existing map. The layer of the rainbows is projected over another layer, that of the clouds, for which we also know the pixel locations. We do this by encoding them as GeoJSON features, which is a kind of description that can be visualised easily on popular mapping technologies like Google maps, d3.js and leaflet.js.
In the app we display the map with the leaflet.js JavaScript library. To get a more austere, black-and-white, look we used the wonderful Open Source map tiles created by Stamen Design called ‘Toner Lite’.
Initially every pixel was encoded as GeoJSON feature. But it turned out this approach was to heavy: the app had to draw thousands of little squares, especially for the clouds. The solution was to group the pixels together in larger shapes: polygons. We accomplished this with Peter Selinger’s potrace, a command line program that can trace bitmap images into vector shapes. Even if potrace has existed for a long time, quit recently Christoph Hormann contributed an output module that allows to create GeoJSON directly. Instead of trying to interpolate smooth curves, we configure potrace as to leave the polygons with straight edges.
