Last week we had the n-th GIS drama about how Mapbox Vector Tiles should be called. I’m actually thiking in creating the GIS version of rubydramas (it’s gone, looks like ruby community has moved to node.js, sorry, io.js). This post is not to talk about naming, we all know standards with company names in the description are never used, look at those ESRI shapefiles…
The objetive of vector tiles (whatever format you use) is to move the data closer to the rendering stage, so it’s projected, clipped, transformed, extra precision is removed, filtered using naive filters, encoded and finally compressed. You should take a look at this Michal Migurski’s blogpost or this talk from Dane Springmeyer in foss4g to understand what is the rationale behind them.
CartoDB is a platform that runs on top of postgis, it renders tiles fetching them directly from a spatial query so we pay the price for all the overhead of going to the database, fetching the data, preparing it for render and so on. That’s why at some point someone though about removing that dynamic part of the equation. But in the other hand we have all the power of postgres and postgis (you know, spatial indices, geometry functions and so on)
So would it possible to generate a easily vector tile from postgres? Of course is, you can do almost everything in postgres right now with an extension but I’m one of those persons that like to do obvious things. It sounds easy, basically we need:
So here it is, given a query that retusn a resultset with cartodb_id and the_geom_webmercator (a geometry column in 3857):
Notice this query does not manage buffer-size, overzooming and so on, that’s pretty easy add tho. Also there is a res/20 that needs an extra explanation. If we used half of the pixel for the snapping we’d soon realize that some polygons and lines are removed pretty soon so using that 20 fixes the thing. I have to say that value was calcualted by hand and there are not maths behind it, why spend hours thinking when with a simple binary search you can fix the thing… The geometry is also simplified after snapping (be sure you do after snapping, the simplify algorithm complexity is higher than the snapping)
does this work?
let’s try with an extract of OSM planet where the geometry is about 350Mb
the CartoDB Vector Tile (did I say I’m pretty good at naming?) is 44Mb (3.8M gzip compressed) so not that bad.
But we still didn’t do anything special with the geometry encoding, we are using WKB to store all the things. Remember that WKB uses 16 bytes per coordinate in a geometry. Mapbox vector tiles use varint encoding of delta values in order to make this smaller. I personally don’t like varint to encode numbers, It’s better to leave the compressor do its work and don’t try to be smart playing with bits. But ok, in postgis we have a way to delta-varint all the things, it’s called TWKB:
copy (select st_astwkbagg(geom, 0, id) from cdb_tile(0, 0, 0, 'select id as cartodb_id, the_geom as the_geom_webmercator from planet')) TO '/tmp/tile2.cvt';
The result is a 9.8Mb (1.8M gzip compressed) tile. Much better and took about the same time to encode it. This also works much better with polygon/lines tables than mixed types, specially when there are a lot of points like in this case.
There a lot of things left, for example, when to use clipping, snapping and simplification (sometimes it’s better to send every single geometry than cut), coincident points, attribute optimization (I didn’t talk about attributes here because with postgres is pretty clear how to do this)
Estaba yo viendo un maravilloso documental sobre los motores turbo en fórmula 1 y me daba cuenta lo parecidas que son todas las áreas de la ingeniería, incluso la informática. El docu dura un par de horas (está formado por dos partes) y básicamente habla del diseño de un motor turbo desde cero con especial énfasis en el uso de la informática como algo novedoso.
Puedes verlo con una palomitas en youtube, incluso si no te interesa para nada la fórmula 1 ni los motores es didáctico y ameno. De hecho me ha asombrado la calidad como documental: riguroso, parándose en los detalles, sin música ni efectos para tener entretenido al personal, sin repetir 200 veces lo mismo antes de hacer una parada de publicidad ni entrevistas a expertos que solo dicen obviedades. En pocas palabras, no tratando como un subnormal al que lo ve.
Y lo mejor, no les ha importado invertir más de media hora de documental enseñándo como trabajaba de verdad esa gente en sus oficinas. No tenéis que buscar mucho para encontrar videos que dan vergüenza ajena sobre como se supone un ingeniero trabaja ahora mismo. Esperad, no busquéis, ya os linko yo uno de twitter
Y es que cuando llegan los retos tecnológicos la realidad del trabajo es como la muestra ese documental. La cagan una y otra vez, Desde el comienzo muestran como van poco avanzando, como se encuentran problemas, como los analizan y solucionan. Joder, incluso piensan.
Me encantan momentos como los siguientes que he vivido unas cuantas veces a lo largo de mi vida profesional:
En el minuto 14 explican como en vez de perder el tiempo diseñando un sistema final usan una solución temporal basada en un compresor externo en vez de un turbo para agilizar las pruebas (os suena?). Ojo al comentario del perico en el minuto 16 cuando les revienta el motor: “we had an accident” (eso sí que es humor del bueno). Más tarde en el 20 se puede ver la cara de desesperación total.
A partir del minuto 48 prueban el nuevo motor que decidieron diseñar desde 0. El motor no arranca, le dan mil vueltas, analizan y al final lo arrancan. Ojo al proceso que siguen y el razonamiento que hacen para solucionarlo (minuto 7 del segundo video).
Cuando ya tiene el motor ya montado en el coche haciendo pruebas en pista, van a probar la parada del motor desde el volante y no funciona (minuto 31, segundo video). Un pequeño detalle estúpido que echa toda la prueba a perder.
El documental está lleno de momentos relacionados con la informática que ahora nos hacen reir pero que en aquel momento eran tecnología punta.
Ahora pienso en como quieren pintar el trabajo que hacemos los que desarrollamos software y me gusta mucho más como es en la realidad, mucho más parecida a como se muestra en ese video. Quiero pensar que cuando este negocio tenga algunos años más nos verán como gente profesional y no como una panda de gilipollas, cortados por el mismo patrón que dan a teclas sin pensar.
CartoDB provides a way to render tiles in non mercator projection, it’s not something we though about when our Map API was designed but it’s possible and works pretty well.
CartoDB stores data in postgres tables using postgis for the geospatial data. I will not explain
this, if you are reading this is because you already know how it works. Tables in CartoDB are
regular tables with some special column names like cartodb_id, the_geom and
the_geom_webmercator
In order to render map tiles CartoDB needs to things:
Maps API uses the_geom_webmercator column by default as geometry source so when a geometry
needs to be rendered you need to call with that name. It also expects that column is projected, no
special SRID is required altough web mercator is for what was designed for.
With that in mind we could tweak Maps API sending a geometry projected in something different than web mercator, so we can create this cartodb.js example:
cartodb.createLayer(map, {
user_name: 'dev',
sublayers: [{
sql: 'select st_transform(the_geom, 32661) as the_geom_webmercator from tm_world_borders_s_11'
cartocss: '#layer { polygon-fill: #F00; polygon-opacity: 0.3; line-color: #F00; }',
}]
})
.addTo(map)
that would render a map projected in SRID 32661, that’s it, it works. The problem is leaflet is still working with mercator. Luckily it provides a way to change the projection used using crs option in the map.
So aplying the same technique than Gregor Aisch describes in this post you can create a CRS that maps the coordinates from the required projection to mercator (and in the other way around) so you can use leaflet as you normally do but using a different projection:
var map = new L.Map('map', {
center: center,
zoom: 2,
crs: cartodb.proj('+proj=stere +lat_0=90 +lat_ts=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs', '32661')
});
cartodb.createLayer(map, {
...
it uses an small library I created for this matter, cartodb.proj
CartoDB not only render png tiles, it also does utf grid tiles to provide interaction and this also works as expected. In the following example an infowindow with some data about the countries is included.
cartodb.createLayer(map, {
user_name: 'dev',
type: 'cartodb',
sublayers: [{
sql: 'select area, iso2, st_transform(the_geom, 32661) as the_geom_webmercator from tm_world_borders_s_11 where st_y(st_centroid(the_geom)) > 0',
cartocss: '#layer { polygon-fill: #F00; polygon-opacity: 0.3; line-color: #F00; }',
interactivity: 'area, iso2'
}]
})
.addTo(map)
.done(function(layer) {
var sub = layer.getSubLayer(0)
cdb.vis.Vis.addInfowindow(map, sub, ['area', 'iso2'])
sub.on('featureOver', function(e, ll, pos, data) {
console.log(data);
})
});
Notice nothing special needs to be added, just cartodb.js code that would work with mercator
Leaflet provides a bunch of methods to work with vector features and guess what, they work as expected so you can use CartoDB SQL API to fetch geometry and render as a GeoJSON layers.
cartodb.SQL({
user: 'dev',
format: 'geojson'
})
.execute('select the_geom from tm_world_borders_s_11 where iso2 = \\'ES\\'')
.done(function(data) {
L.geoJson(data).addTo(map);
});
L.marker(center).addTo(map); </code>
in this case the projection is done client side, notice the example fetchs 4326 from CartoDB
Some projections don’t work well depending on the zone you are showing and that’s something you should expect, most projections are meant to work well only in a concrete region. For example, a projection ready to show Antartica is not going to work well showing Norway and it will be rendered with distorsion or it will not even rendered because the projection fails.
In PostGIS some projections fails when you try to render outside a bounding box:
db=# select st_transform(st_setsrid(st_makepoint(0, -181), 4326), 2857);
ERROR: transform: couldn't project point (0 -181 0): latitude or longitude exceeded limits (-14)
in order to make it work we can use PostGIS functions ST_Intersection, like:
SELECT ST_Intersection(the_geom, ST_MakeEnvelope(-180, 0, 180, 90, 4326))
in this case we get the geometry in the north side or the earth.
A map I created with CartoDB with some dataset Antartica
Faltan aún unos días para cerrar el año pero voy a adelantar el post para ver si lo cerramos un poco antes. No me extenderé demasiado este año.
En mi mesilla de noche tengo una foto que me hice el día que presenté mi proyecto fin de carrera con mis abuelos. Es lo último que veo nítido antes de quitarme las gafas de hipster que me compré hace unos meses.
La foto es bastante mala pero me hace recordar dos cosas: a mi abuela que se fue este año sin avisar y los valores que aprendí de ellos.
Este año he hecho el gilipollas a nivel que ni yo mismo me creo, he pasado las líneas que marcan esos valores con mucho, así que cuando me voy a la cama ese “hasta mañana” a la foto me ayudar a recordar donde están los esos límites que me enseñaron.
One of the things I’m doing lately is analyze datasets in order to improve rendering speed in CartoDB and see how postgres performs. Unfortunately (or not, who knows) I’m don’t know so much about geospatial indices so I try to analyze using a bunch of SQL queries with the help of explain analyze.
On our tiler we have a fixed number of mapnik workers that can run at the same time (obviously) and setting up a new worker to render an empty tile is too expensive or at least more than executing a SQL query to see if there is data in this tile.
I was thinking about how would be to generate a index to khow, given a SQL query, what are the empty tiles. One of the thinks that came to my mind was the recursive WITH postgres statement. At the end webmercator is based on a quadtree so we can iterate it using recursion. This is the query I wrote:
WITH RECURSIVE t(x, y, z, e) AS (
-- root node (0, 0, 0)
SELECT 0, 0, 0, exists(select 1 from ships where the_geom_webmercator && CDB_XYZ_Extent(0, 0, 0))
UNION ALL
-- coordinate for the children
SELECT x*2 + xx, y*2 + yy, z+1,
exists(select 1 from ships where the_geom_webmercator && CDB_XYZ_Extent(x*2 + xx, y*2 + yy, z+1)) from t,
-- iterate over 4 children
(VALUES (0, 0), (0, 1), (1, 1), (1, 0)) as c(xx, yy)
-- only for tiles with geometry and up to zoom level 8
where e AND z < 8
)
SELECT z, x, y FROM t where e
(CDB_XYZ_Extent method returns the bbox for the tile (x, y, z))
It takes 2.8 seconds on my laptop for a dataset distributed all over the world with ~800k points
This generates a table where you can lookup if a tile should be rendered.
Another similar approach is to analyze what tiles you should fetch in order to not retrieve too much information (lets say more than 4k points) but not query the database for tiles with a few points when you can use some parent tile to render (and use cache for that). The SQL query is similar, the only thing that changes is the count(*)
WITH RECURSIVE t(x, y, z, e) AS (
SELECT 0, 0, 0, (select count(*) from ships where the_geom_webmercator && CDB_XYZ_Extent(0, 0, 0))
UNION ALL
SELECT x*2 + xx, y*2 + yy, z+1, (select count(*) from ships where the_geom_webmercator && CDB_XYZ_Extent(x*2 + xx, y*2 + yy, z+1)) from t, (VALUES (0, 0), (0, 1), (1, 1), (1, 0)) as c(xx, yy) where e > 4000 AND z < 17
)
SELECT z, x, y, e FROM t where e > 0
notice in this case I raised the limit to 17, the recursion in the case stops before due the number of points, there are still 37 tiles with more than 4k points tho. It takes 9.8 seconds to generate for the same dataset. If we increase the number of points to 64k the time is reduced to a half.
I also created a map to see the quadtree in action: