Tuesday, February 25, 2014

3D terrain visualization with python and Mayavi2

I have always wanted to draw these 3D terrains like those in www.shadedrelief.com, which are amazing. But the examples were all using software I don't use, so I tried to do it with python.
 The final result

As usual, you can get all the source code and data at my GitHub page.

Getting the data

After trying different locations, I decided to use the mountain of Montserrat, close to Barcelona, since it has nice stone towers that are a good test for the DEM and the 3D visualization. An actual picture of the zone used is this one:
 Montserrat monastery
The building is a good reference, since the stone only areas make the result testing much more difficult.
All the data has been downloaded from the ICGC servers:
• The DEM data was downloaded from the Vissir3 service, going to the section catàleg i descàrregues -> MDE 5x5. The file is named met5v10as0f0392Amr1r020.txt, but I cut a small part of it, to make mayavi2 work smoother using:

gdalwarp -te 401620 4604246 403462 4605867 -s_srs EPSG:25831 -t_srs EPSG:25831 met5v10as0f0392Amr1r020.txt dem.tiff
• The picture to drap over the dem file can be downloaded using the WMS service given by the ICGC:

http://geoserveis.icc.cat/icc_mapesbase/wms/service?REQUEST=GetMap&VERSION=1.1.0&SERVICE=WMS&SRS=EPSG:25831&BBOX=401620.0,4604246.0,403462.0,4605867.0&WIDTH=1403&HEIGHT=1146&LAYERS=orto5m&STYLES=&FORMAT=JPEG&BGCOLOR=0xFFFFFF&TRANSPARENT=TRUE&EXCEPTION=INIMAGE
It's not as automatic as I would like, but if it's possible to download a DEM and the corresponding image, it's possible to create the 3D image.

Creating the image

First, let's plot the DEM file in 3D using mayavi2:
"""
Plotting the terrain DEM with Mayavi2
"""

from osgeo import gdal
from mayavi import mlab

ds = gdal.Open('dem.tiff')

mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))

mlab.surf(data, warp_scale=0.2)

mlab.show()
• In first place, we import as usual, gdal and numpy. Also the mlab library from mayavi, which lets set the mayavi canvas.
• The data is read, as usual, with the gdal ReadAsArray method.
• The figure is created. This works like creating the Image object in the PIL library, creating the canvas where the data wil be drawn. In this case, the size is 640 x 800 pixels, making the figure bigger can affect the performance in some old computers. bgcolor sets the blue color as the background.
• The surf method will plot the surface. The input has to be a 2D numpy array, which is what we have.
• The warp_scale argument sets the vertical scale. In this case, letting the default value (1?) creates a really exaggerated effect so its better to play a little to get a more realistic effect.
• The colors depend of the Z value at each point, and can be changed using the color or colormap option.
• The show() method makes the image to stay visible when running the example from a script. If you use ipython, you don't need this step.
• If you want to save the figure as a png, you can either use the icon in the mayavi window or call the method mlab.savefig('image_name')
• If you want to move the camera (change the prespective), you can use the roll/yaw/pitch methods:
f = mlab.gcf()
camera = f.scene.camera
camera.yaw(45)

 The plotted DEM
Now, let's put an aerial image over the 3D visualization:
"""
Draping an image over a terrain surface
"""
from osgeo import gdal
from tvtk.api import tvtk
from mayavi import mlab
import Image

ds = gdal.Open('dem.tiff')
im1 = Image.open("ortofoto.jpg")
im2 = im1.rotate(90)
im2.save("/tmp/ortofoto90.jpg")
bmp1.file_name="/tmp/ortofoto90.jpg" #any jpeg file

my_texture=tvtk.Texture()
my_texture.interpolate=0
my_texture.set_input(0,bmp1.get_output())

mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))

surf = mlab.surf(data, color=(1,1,1), warp_scale=0.2)
surf.actor.enable_texture = True
surf.actor.tcoord_generator_mode = 'plane'
surf.actor.actor.texture = my_texture

mlab.show()

• The most important new import is tvtk. TVTK is a python api that allows to work with VTK objects. Actually, my knowledge of Mayavi2 is very limited, but I see TVTK as an extension.
• The DEM data is read the same way, using the ReadAsArray method.
• The aerial image, named ortofoto.jpg, is not in the correct orientation. It took me a lot of time to get what was happening. I rotate it opening the image with the PIL library and using the rotate method (lines 11 to 13)
• Then, the tvtk object with the texture is created, loading the image with a JPEGReader object, and assigning it to the Texture object (lines 14 to 19).
• The figure and the 3D surface is created as in the other example (lines 22 and 24)
• Then, the surface is modified to show the image over it (lines 25 to 28).
 The final result

The result is correct, but the aerial image is a little moved from the place it should be so, since the terrain is really steep, part of the buildings are drawn on the stone walls! I edited the WMS coordinates a little so the result is slightly better. Anyway, the method is correct.

• GitHub: The source code and example data
• Mayavi2: 3D scientific data visualization and plotting
• This entry to a mailing list gave me the tip to create the example
• If you want to do more or less the same, but using JavaScript, Bjørn Sandvik posted this excellent example.
• ShadedRelief: Cool 3D and shaded relief examples.

Saturday, January 18, 2014

D3 map Styling tutorial III: Drawing animated paths

The example La Belle France, or the original La Bella Italia by Gregor Aisch, have a nice ferry sailing along a path to connect the islands with the continent.
Drawing a path in a map, and animating some icon on it can be a nice tool to show information about routes, storm tracks, and other dynamic situations.

This example shows how to draw the Haiyan typhoon track on the map drawn in the last post.

The working examples are here:
Creating a geopath
Animating an object on the path

Getting the data

Both the base map data and the typhoon data are explained in the post  D3 map Styling tutorial I: Preparing the data

Creating a geopath

First, how to draw the path line on a map. The working example is here.



.graticule {
fill: none;
stroke: #777;
stroke-opacity: .5;
stroke-width: .5px;
}

.land {
fill: #999;
}

.boundary {
fill: none;
stroke: #fff;
stroke-width: .5px;
}

svg .path {
fill: none;
stroke-opacity: .8;
stroke-dasharray: 3,2;
stroke: #f44;
}

var width = 600,
height = 500;

var projection = d3.geo.mercator()
.scale(5*(width + 1) / 2 / Math.PI)
.translate([width / 2, height / 2])
.rotate([-125, -15, 0])
.precision(.1);

var path = d3.geo.path()
.projection(projection);

var graticule = d3.geo.graticule();

var svg = d3.select("body").append("svg")
.attr("width", width)
.attr("height", height);

svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("d", path);
d3.json("track.json", function(error, track) {
d3.json("/mbostock/raw/4090846/world-50m.json", function(error, world) {
svg.insert("path", ".graticule")
.datum(topojson.feature(world, world.objects.land))
.attr("class", "land")
.attr("d", path);

svg.insert("path", ".graticule")
.datum(topojson.mesh(world, world.objects.countries, function(a, b) { return a !== b; }))
.attr("class", "boundary")
.attr("d", path);

var pathLine = d3.svg.line()
.interpolate("cardinal")
.x(function(d) { return projection([d.lon, d.lat])[0]; })
.y(function(d) { return projection([d.lon, d.lat])[1]; });

var haiyanPath = svg.append("path")
.attr("d",pathLine(track))
.attr("class","path");

});
});

d3.select(self.frameElement).style("height", height + "px");


• The base map is drawn in the simplest way, as shown in this example, so the script stays clearer.
• The typhoon track is loaded from the json file generated in the first tutorial post (line 55)
• The path is created and inserted from lines 68 to 75:
• A d3.svg.line element is created. This will interpolate a line between the points. An other option is to draw segments from each point, so the line is not so smooth, but the actual points are more visible.
• The interpolate method sets the interpolation type to be used.
• x and y methods, set the svg coordinates to be used. In our case, we will transform the geographical coordinates using the same projection function set for the map. The coordinates transformation is done twice, one for the x and another for the y. It would be nice to do it only once.
• The path is added to the map, using the created d3.svg.line, passing the track object as a parameter to be used by the line function. The class is set to path, so is set to a dashed red line (line 20)
Drawing the paths is quite easy, taking only two steps.

Animating an object on the path

The typhoon position for every day is shown on the path, with an icon. The icon size and color change with the typhoon class. The working example is here.


.graticule {
fill: none;
stroke: #777;
stroke-opacity: .5;
stroke-width: .5px;
}

.land {
fill: #999;
}

.boundary {
fill: none;
stroke: #fff;
stroke-width: .5px;
}

var width = 600,
height = 500;

var projection = d3.geo.mercator()
.scale(5*(width + 1) / 2 / Math.PI)
.translate([width / 2, height / 2])
.rotate([-125, -15, 0])
.precision(.1);

var path = d3.geo.path()
.projection(projection);

var graticule = d3.geo.graticule();

var svg = d3.select("body").append("svg")
.attr("width", width)
.attr("height", height);

svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("d", path);
d3.json("track.json", function(error, track) {
d3.json("/mbostock/raw/4090846/world-50m.json", function(error, world) {
var color_scale = d3.scale.quantile().domain([1, 5]).range(colorbrewer.YlOrRd[5]);

var filter = svg.append("defs")
.append("filter")
.attr("height", "130%");
filter.append("feGaussianBlur")
.attr("in", "SourceAlpha")
.attr("stdDeviation", 5)
.attr("result", "blur");

filter.append("feOffset")
.attr("in", "blur")
.attr("dx", 5)
.attr("dy", 5)
.attr("result", "offsetBlur");

var feMerge = filter.append("feMerge");

feMerge.append("feMergeNode")
.attr("in", "offsetBlur")
feMerge.append("feMergeNode")
.attr("in", "SourceGraphic");

svg.insert("path", ".graticule")
.datum(topojson.feature(world, world.objects.land))
.attr("class", "land")
.attr("d", path)

svg.insert("path", ".graticule")
.datum(topojson.mesh(world, world.objects.countries, function(a, b) { return a !== b; }))
.attr("class", "boundary")
.attr("d", path);

var dateText = svg.append("text")
.attr("id", "dataTitle")
.text("2013/11/"+track[0].day + " " + track[0].hour + ":00 class: " + track[0].class)
.attr("x", 70)
.attr("y", 20)
.attr("font-family", "sans-serif")
.attr("font-size", "20px")
.attr("fill", color_scale(track[0].class));

var pathLine = d3.svg.line()
.interpolate("cardinal")
.x(function(d) { return projection([d.lon, d.lat])[0]; })
.y(function(d) { return projection([d.lon, d.lat])[1]; });

var haiyanPath = svg.append("path")
.attr("d",pathLine(track))
.attr("fill","none")
.attr("stroke", color_scale(track[0].class))
.attr("stroke-width", 3)

.style('stroke-dasharray', function(d) {
var l = d3.select(this).node().getTotalLength();
return l + 'px, ' + l + 'px';
})
.style('stroke-dashoffset', function(d) {
return d3.select(this).node().getTotalLength() + 'px';
});

var haiyanPathEl = haiyanPath.node();
var haiyanPathElLen = haiyanPathEl.getTotalLength();

var pt = haiyanPathEl.getPointAtLength(0);

var icon = svg.append("path")
.attr("d","m 20,-42 c -21.61358,0.19629 -34.308391,10.76213 -41.46346,18.0657 -7.155097,7.3036 -11.451337,17.59059 -11.599112,26.13277 0,14.45439 9.037059,26.79801 21.767213,31.69368 -14.965519,10.64929 -25.578236,6.78076 -37.671451,7.85549 C -4.429787,54.20699 14.03,37.263 23.12144,28.41572 32.2133,19.56854 34.6802,10.79063 34.82941,2.19847 c 0,-14.45219 -9.03405,-26.79679 -21.76113,-31.69364 14.90401,-10.54656 25.48889,-6.69889 37.55061,-7.77104 C 38.78869,-40.57565 29.11666,-41.95733 21.03853,-42 20.68954,-42.0105 20.34303,-42.0105 20,-42 z M 0.82306,-7.46851 c 4.72694,0 8.56186,4.27392 8.56186,9.54602 0,5.2725 -3.83492,9.54651 -8.56186,9.54651 -4.726719,0 -8.555958,-4.27401 -8.555958,-9.54651 0,-5.2721 3.829239,-9.54602 8.555958,-9.54602 z")
.attr("transform", "translate(" + pt.x + "," + pt.y + "), scale("+(0.15*track[0].class)+")")
.attr("fill", color_scale(track[0].class))
.attr("class","icon");

var i = 0;
var animation = setInterval(function(){
pt = haiyanPathEl.getPointAtLength(haiyanPathElLen*i/track.length);
icon
.transition()
.ease("linear")
.duration(1000)
.attr("transform", "translate(" + pt.x + "," + pt.y + "), scale("+(0.15*track[i].class)+"), rotate("+(i*15)+")")
.attr("fill", color_scale(track[i].class));

haiyanPath
.transition()
.duration(1000)
.ease("linear")
.attr("stroke", color_scale(track[i].class))
.style('stroke-dashoffset', function(d) {
var stroke_offset = (haiyanPathElLen - haiyanPathElLen*i/track.length + 9);
return (haiyanPathElLen < stroke_offset) ? haiyanPathElLen : stroke_offset + 'px';
});

dateText
.text("2013/11/"+track[i].day + " " + track[i].hour + ":00 class: " + track[i].class)
.attr("fill", color_scale(track[i].class));
i = i + 1;
if (i==track.length)
clearInterval(animation)

},1000);

});
});

d3.select(self.frameElement).style("height", height + "px");


This second example is more complex than the first one:
• The base map has a shadow effect. See the second part of the tutorial for the source.
• The map is animated:
• Line 135 sets an interval, so the icon and line can change with the date.
• A variable i is set, so the array elements are used in every interval.
• When the dates have ended, the interval is removed, so everything stays quiet. Line 158.
• An icon moves along the path indicating the position of the typhoon
•  Line 128 created the icon. First, I created it using inkscape, and with the xml editor that comes with it, I copied the path definition. This method can get really complex with bigger shapes.
• Line 136 finds the position of the typhoon. The length of the track is found at line 122 with the getTotalLength() method.
• Line 137 moves the icon. A transition is set, so the movement is continuous even thought the points are separated. The duration is the same as the interval, so when the icon has arrived at the final point, a new transformation starts to the next one.
• Line 141 has the transform operation that sets the position (translate), the size (scale) and rotation (rotate). The factors multiplying the scale and rotation are those only to adjust the size and rotation speed. They are completely arbitrary.
• The path gets filled when the icon has passed. I made this example to learn how to do it. Everything happens at line 144. Basically, the trick is creating a dashed line, and playing with the stroke-dashoffset attribute to set where the path has to arrive.
• The color of the path and icon change with the typhoon class
• At line 54, a color scale is created using the method d3.scale.quantile
• The colors are chosen with colorbrewer, which is a set of color scales for mapping, and has a handy javascript library to set the color scales just by choosing their name. I learned how to use it with this example by Mike Bostock.
• The lines 156 and 142 change the track and icon colours.
• Finally, at line 154, the date is changed, with the same color as the typhoon and track.

Simple path on a map - The first example
Haiyan typhoon track - The second example
D3 map Styling tutorial I: Preparing the data
D3 map Styling tutorial II: Giving style to the base map
Animated arabic kufic calligraphy with D3 - Animating paths using d3js
La Bella Italia - Kartograph example by Gregor Aisch
La Belle France - The same example as La Bella Italia, but using D3js
Every ColorBrewer Scale - An example to learn how to use ColorBrewer with D3js

Saturday, January 4, 2014

D3 map Styling tutorial II: Giving style to the base map

The example La Belle France, or the original La Bella Italia by Gregor Aisch use svg filters to give style to the maps. I haven't found many examples about how to do it, so I will explain here two different styles: a simple shadow under the map, and a style similar to the original map.

As usual, the examples can be seen at my bl.ocks.org page:
LaBella Italia

SVG filters basics

SVG can style the elements by using css or css-like attributes, like the color, the stroke, and so on. But it's possible to use filters to add efects such as blurring, dilating the shapes, adding lights... This, of course, can be used to change the map polygons too.

A basic usage of an SVG filter is like this:



And the result is:
Note that:

• The shape (a rectangle defined by the rect tag) is drawn as usual
• The filter is defined with the filter tag inside the defs section. The filter must have an id to be used in the geometries where it has to be used
• The filter is applied to the rectangle by using the filter="url(#blur)" attribute, where the id of the filter is added as url(#name_of_the_filter)
This is the easiest filter an SVG can have, but filters can concatenate different effects one after the other or in parallel, as we will see in the following examples.
A list of all the effects and a tutorial can be found at W3C Schools

The drop shadow effect gives a quite nice look to the maps and it's fast to process by the browser and easy to code. I took the code from this example.

The map shown in the last post with the effect. The code:


.graticule {
fill: none;
stroke: #777;
stroke-opacity: .5;
stroke-width: .5px;
}

.land {
fill: #999;
}

.boundary {
fill: none;
stroke: #fff;
stroke-width: .5px;
}

var width = 600,
height = 500;

var projection = d3.geo.mercator()
.scale(5*(width + 1) / 2 / Math.PI)
.translate([width / 2, height / 2])
.rotate([-125, -15, 0])
.precision(.1);

var path = d3.geo.path()
.projection(projection);

var graticule = d3.geo.graticule();

var svg = d3.select("body").append("svg")
.attr("width", width)
.attr("height", height);

var filter = svg.append("defs")
.append("filter")
.attr("height", "130%");
filter.append("feGaussianBlur")
.attr("in", "SourceAlpha")
.attr("stdDeviation", 5)
.attr("result", "blur");

filter.append("feOffset")
.attr("in", "blur")
.attr("dx", 5)
.attr("dy", 5)
.attr("result", "offsetBlur");

var feMerge = filter.append("feMerge");

feMerge.append("feMergeNode")
.attr("in", "offsetBlur")
feMerge.append("feMergeNode")
.attr("in", "SourceGraphic");

svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("d", path);
d3.json("/mbostock/raw/4090846/world-50m.json", function(error, world) {
svg.insert("path", ".graticule")
.datum(topojson.feature(world, world.objects.land))
.attr("class", "land")
.attr("d", path)

svg.insert("path", ".graticule")
.datum(topojson.mesh(world, world.objects.countries, function(a, b) { return a !== b; }))
.attr("class", "boundary")
.attr("d", path);

});

d3.select(self.frameElement).style("height", height + "px");


Note that the filter (lines 46 to 66) has several parts:

• The filter is defined adding a defs tag and appending the filter, with an id and a height tags. The height is added to affect the region outside the geometry. It's not actually necessary in this case, but I kept it to maintain the example code.
• An feGaussianBlur is then appended (line 50). The in attribute indicates that the input image for the filter is the alpha value of the input image, so if the map is colored in red, the shadow will remain grey. stdDeviation indicates the filter intensity, and output gives an id to use the resulting image, as we will see.
• Line 55 adds an feOffset filter, that will move the shadow. The dx and dy attributes define how to move the shadow, but the important attribute here is the in, that takes the id defined in the feGaussianBlur tag, so what is moved is the blurred image. Again, an out attribute is defined to use the result image later.
• Line 61 defines an feMerge tag. This is an interesting property of svg filters that allow to append different filter results one over the other to create more complex outputs. Appended to the tag, an feMergeNode is added at line 63 with the blurred image at the in attribute. At line 65, the original image is appended over the blurred image, with an other feMergeNode tag, with the in attr set to SourceGraphic to indicate the original image.
SVG filters seemed quite dificult to me, but they are actually a kind of instructions one after the other, written in XML.

Styling like La Bella Italia

As it was shown in the post, the filters part of La Bella Italia map is quite complicated. I'll try to show an easy way to get a similar effect.
Disclaimer: In my old computer, generating these maps is very slow, so maybe it would be better to think about another combination to create cool maps without using the erode and dilate filters.

The map without any effect looks like this:
And the effect has three stacked parts (I'll paste only the part of the code affecting each part, the whole code is here):

1. A colouring and blurring
To do it:
filter.append("feColorMatrix")
.attr("in","SourceGraphic")
.attr("type", "matrix")
.attr("values", "0 0 0 0 0.6 0 0 0 0 0.5333333333333333 0 0 0 0 0.5333333333333333  0 0 0 1 0")
filter.append("feGaussianBlur")
.attr("stdDeviation", 15)
.attr("result", "f1blur");
So there are two concatenated actions. The first one, feColorMatrix, changes the color of the original image (see how here). The second, blurs it as in the first example.
2. Eroding, coloring, blurring and composing:
filter.append("feColorMatrix")
.attr("in","SourceGraphic")
.attr("type", "matrix")
.attr("values", "0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 500 0")
filter.append("feMorphology")
.attr("operator","erode")
.attr("result","f2r1");
filter.append("feGaussianBlur")
.attr("in","f2r1")
.attr("stdDeviation","4")
.attr("result","f2r2");
filter.append("feColorMatrix")
.attr("in","f2r2")
.attr("type", "matrix")
.attr("values", "1 0 0 0 0.5803921568627451 0 1 0 0 0.3607843137254902 0 0 1 0 0.10588235294117647 0 0 0 -1 1")
.attr("result","f2r3");
filter.append("feComposite")
.attr("operator","in")
.attr("in","f2r3")
.attr("result","f2comp");
This one is more complicated.
1. The first step changes the map color into black, with an alpha value of 0.5. This output will be used in the second and last step.
2. Then, using the first output, the image is eroded. That is, the land gets smaller by one pixel, using feMorphology with the erode operator. Then, the result is blurred using feGaussianBlur, and coloured as in the first filter.
3. The resulting image is composited with the original black and white image, using feComposite. The definition of this operator, taken from here, is:
The result is the part of A that is within the boundaries of B. Don't confuse the name of this attribute value with the in attribute.
3. The two effects are stacked, and the original map is added at the end:
var feMerge = filter.append("feMerge");

feMerge.append("feMergeNode")
.attr("in", "f1blur");
feMerge.append("feMergeNode")
.attr("in", "f2comp");
feMerge.append("feMergeNode")
.attr("in", "SourceGraphic");
This part is simple,the feMerge stacks all the outputs indicated in the feMergeNode tags. So in this case, the blurred image goes in the first place, then the eroded one, and finally, the original one. Now, the effect seems quite a lot to the one in La Bella Italia.

Drop shadow map - The first example
La Bella Italia like map - The second example
D3 map Styling tutorial I: Preparing the data - First part of the tutorial
La Bella Italia - Kartograph example by Gregor Aisch
La Belle France - The same example as La Bella Italia, but using D3js
d3.js drop shadow example - SVG drop shadow filter by Charl P. Botha
SVG filters tutorial - W3C Schools tutorial about SVG filters

Monday, December 30, 2013

D3 map Styling tutorial I: Preparing the data

This is the first post of a series to explain how to draw styled maps with paths in them, using D3js.
The final map will be an animated map of the Haiyan typhoon track:
This entries are to explain my last post La Belle France, which had a part with filters and an other with paths, not very clearly explained.
So first, the typhoon path data is needed, and a base map to draw under it too.

Base Map

To draw the maps for this tutorial, the wolrd map at 50m resolution used in most of the mbostock's examples will be used. The base map, is taken from this example, and can be seen, as usual, at bl.ocks.org.


.graticule {
fill: none;
stroke: #777;
stroke-opacity: .5;
stroke-width: .5px;
}

.land {
fill: #999;
}

.boundary {
fill: none;
stroke: #fff;
stroke-width: .5px;
}

var width = 600,
height = 500;

var projection = d3.geo.mercator()
.scale(5*(width + 1) / 2 / Math.PI)
.translate([width / 2, height / 2])
.rotate([-125, -15, 0])
.precision(.1);

var path = d3.geo.path()
.projection(projection);

var graticule = d3.geo.graticule();

var svg = d3.select("body").append("svg")
.attr("width", width)
.attr("height", height);

svg.append("path")
.datum(graticule)
.attr("class", "graticule")
.attr("d", path);
d3.json("/mbostock/raw/4090846/world-50m.json", function(error, world) {
svg.insert("path", ".graticule")
.datum(topojson.feature(world, world.objects.land))
.attr("class", "land")
.attr("d", path);

svg.insert("path", ".graticule")
.datum(topojson.mesh(world, world.objects.countries, function(a, b) { return a !== b; }))
.attr("class", "boundary")
.attr("d", path);

});

d3.select(self.frameElement).style("height", height + "px");



The basic difference from the Mercator example is the projection variable:
var projection = d3.geo.mercator()
.scale(5*(width + 1) / 2 / Math.PI)
.translate([width / 2, height / 2])
.rotate([-125, -15, 0])
.precision(.1);

• The scale is multoplied by five, to zoom to the typhoon zone
• A rotation is added, to center the map to the typhoon zone

Typhoon Haiyan track

It's not easy to find the actual data for the hurricanes and typhoons in a convenient format. After looking into several sites, the one I like the most is this one from Japan. The data is in a list which is possible to copy, but putting all the points into an array is not an amusing activity. So here is the script I've used to transform the data into a json:
f = open("track.txt","r")
out = "["
for line in f:
day = int(line[12:16])
hour = int(line[16:18])
lat = float(line[20:24])
lon = float(line[28:34])
h_class = int(line[48:49])

out += '{"day":%d, "hour":%d, "lat":%.1f, "lon":%.1f, "class": %d},' %(
day, hour, lat, lon, h_class)
out = out[0:-1] + "]"
f.close()

print out
Where track.txt is the file where I pasted the data copied from the web page.

The base map for the next posts
La Belle France - The blog entry that made me write this post
D3 Mercator world map example
Haiyan typhoon track data, from the Japan Meteorological Agency Best Track Data
Digital Typhoon in English, with tracks from many cuurent and past typhoons
Hong Kong observatory Haiyan typhoon track data

Tuesday, December 3, 2013

La Belle France - Map styling with D3js

When I looked to the vectorial web mapping tools for the first time, I was completely carried away by this example made with Kartograph by the software creator Gregor Aisch.
This example was, actually, the reason I learned Kartograph and published some posts here.
Later, I learned D3js, which was also amazing, but I didn't find a map like La Bella Italia. So I tried to do a similar one myself.

To make the map, I had to reproduce a the result of the Kartograph example to understand how the effects were achieved. Then, get some base data to make the map, and finally, code the example.

Getting the data

The map has, basically, data about the land masses, the countries and the French regions.
I tried data from different places (one of them, Eurostat, that ended in this example), until I decided to use the Natural Earth data. After some attempts, I decided to use the 1:10m data. The files are:

1. ne_10m_admin_0_countries.shp, clipped using:
ogr2ogr -clipsrc -10 35 15 60 countries.shp ne_10m_admin_0_countries.shp
2. ne_10m_admin_1_states_provinces.shp, clipped using:
ogr2ogr -clipsrc -10 35 15 60 -where "adm0_a3 = 'FRA'" regions.shp ne_10m_admin_1_states_provinces.shp
3. ne_10m_land.shp, downloaded from github, since the official version gave errors when converted to TopoJSON. Clipped using:
ogr2ogr -clipsrc -10 35 15 60  land.shp ne_10m_land.shp
With that, the land, countries and regions are available. To merge them into a single TopoJson file, I used:
topojson -o ../data.json countries.shp regions.shp

The html code

Since the code is quite long, and I think that I will made some more posts about specific parts of the technique, the comments are a bit shorter than usually.

CSS

Note how is the font AquilineTwo.ttf  loaded:
@font-face {
font-family: 'AquilineTwoRegular';
src: url('AquilineTwo-webfont.eot');
src: url('AquilineTwo-webfont.eot?#iefix') format('embedded-opentype'),
url('AquilineTwo-webfont.woff') format('woff'),
url('AquilineTwo-webfont.ttf') format('truetype'),
url('AquilineTwo-webfont.svg#AquilineTwoRegular') format('svg');
font-weight: normal;
font-style: normal;

}


Later, the font can be set using .attr("font-family","AquilineTwoRegular")

To achieve the effects, some layers are loaded more than one time, so different filters can be applied to get the shades and blurred borders:
svg.selectAll(".bgback")
.data(topojson.feature(data, data.objects.land).features)
.enter()
.append("path")
.attr("class", "bgback")
.attr("d", path)
.style("filter","url(#oglow)")
.style("stroke", "#999")
.style("stroke-width", 0.2);
In this case, the land masses are drawn, applying the effect named oglow, which looks like:
var oglow = defs.append("filter")
.attr("id","oglow");
oglow.append("feColorMatrix")
.attr("in","SourceGraphic")
.attr("type", "matrix")
.attr("values", "0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 1 0")
oglow.append("feMorphology")
.attr("operator","dilate")
oglow.append("feColorMatrix")
.attr("type", "matrix")
.attr("values", "0 0 0 0 0.6 0 0 0 0 0.5333333333333333 0 0 0 0 0.5333333333333333  0 0 0 1 0")
.attr("result","r0");
oglow.append("feGaussianBlur")
.attr("in","r0")
.attr("stdDeviation","4")
.attr("result","r1");
oglow.append("feComposite")
.attr("operator","out")
.attr("in","r1")
.attr("result","comp");
To see how svg filters work, many pages are available. I got them looking at the Kartograph example generated html.

The labels aren't inside the TopoJSON (although they could be!), so I decided the labels to add and put them into an array:
var cities = [
{'pos': [2.351, 48.857], 'name': 'Paris'},
{'pos':[5.381, 43.293], 'name': 'Marseille'},
{'pos':[3.878, 43.609], 'name': 'Montpellier'},
{'pos':[4.856, 45.756], 'name': 'Lyon'},
{'pos':[1.436, 43.602], 'name': 'Toulouse'},
{'pos':[-0.566, 44.841], 'name': 'Bordeaux'},
{'pos':[-1.553, 47.212], 'name': 'Nantes'},
{'pos':[8.737, 41.925], 'name': 'Ajaccio'},
];
Then, adding them to the map is easy:
var city_labels =svg.selectAll(".city_label")
.data(cities)
.enter();

city_labels
.append("text")
.attr("class", "city_label")
.text(function(d){return d.name;})
.attr("font-family", "AquilineTwoRegular")
.attr("font-size", "18px")
.attr("fill", "#544")
.attr("x",function(d){return projection(d.pos)[0];})
.attr("y",function(d){return projection(d.pos)[1];});

city_labels
.append("circle")
.attr("r", 3)
.attr("fill", "black")
.attr("cx",function(d){return projection(d.pos)[0];})
.attr("cy",function(d){return projection(d.pos)[1];});
Note that the positions must be calculated transforming the longitude and latitude using the d3js projection functions.

The ship

To draw the ship, tow things are necessary, the path and the ship.
To draw the path:
var ferry_path = [[8.745, 41.908],
[8.308, 41.453],
[5.559, 43.043],
[5.268, 43.187],
[5.306, 43.289]
];
var shipPathLine = d3.svg.line()
.interpolate("cardinal")
.x(function(d) { return projection(d)[0]; })
.y(function(d) { return projection(d)[1]; });

var shipPath = svg.append("path")
.attr("d",shipPathLine(ferry_path))
.attr("stroke","#000")
.attr("class","ferry_path");
Basically, d3.svg.line is used to interpolate the points, making the line smoother. This is easier than the Kartograph way with geopaths, where the Bézier control points have to be calculated. d3.svg.line is amazing, more than what I thought before.
I don't know if the way to calculate the projected points is the best one, since I do it twice for each point, which is ugly.
To move the ship, a ship image is appended, and then moved with a setInterval:
  var shipPathEl = shipPath.node();
var shipPathElLen = shipPathEl.getTotalLength();

var pt = shipPathEl.getPointAtLength(0);
var shipIcon = svg.append("image")
.attr("x", pt.x - 10)
.attr("y", pt.y - 5.5)
.attr("width", 15)
.attr("height", 8);

var i = 0;
var delta = 0.05;
var dist_ease = 0.2;
var delta_ease = 0.9;
setInterval(function(){

pt = shipPathEl.getPointAtLength(i*shipPathElLen);
shipIcon
.transition()
.ease("linear")
.duration(1000)
.attr("x", pt.x - 10)
.attr("y", pt.y - 5.5);

//i = i + delta;

if (i < dist_ease){
i = i + delta * ((1-delta_ease) + i*delta_ease/dist_ease);
}else if (i > 1 - dist_ease){
i = i + delta * (1 - ((i - (1 - dist_ease)) * (delta_ease/dist_ease)));
}else{
i = i + delta;
}
if (i+0.0001 >= 1 || i-0.0001 <= 0)
delta = -1 * delta;
},1000);
The ship position is calculated every second, and moved with a d3js transition to make it smooth (calculating everything more often didn't give this smooth effect)
The speed of the ship is changed depending on the proximity to the harbour, to a void the strange effect of the ship crashing into it. The effect is controlled by dist_ease and delta_ease parameters, that change the distance where the speed is changed, and the amount of speed changed.

What's next

• The SVG filters should be explained in  a better way, maybe packing them into functions as Kartograph does.
• SVG rendering lasts quite a lot in my computer. The same happens with Kartograph, so the problem comes from the SVG rendering. Anyway, could be improved.
• A canvas version would be nice.