[OSM-talk] mapnik db rendering

[Artem Pavlenko]

On Wednesday 17 January 2007 15:36, SteveC wrote:
> actually no, I'm going to give up as I don't know what colours things
> are supposed to be.
>
> Richard, or anyone, what way round should the pri/sec/trunk roads be?
>

OK, I've updated osm.xml in SVN trunk. This is how it looks now:
http://media.mapnik.org/images/enfield.jpg
Please, can anyone let me know if this is correct!

> Artem, how do you convert from google zoom <--> mapnik scale?
Assuming you're after Mapnik's scale demoninators corresponding to GYM 
levels , here it is:
zoom level=1 scale_denom=279541132.014
zoom level=2 scale_denom=139770566.007
zoom level=3 scale_denom=69885283.0036
zoom level=4 scale_denom=34942641.5018
zoom level=5 scale_denom=17471320.7509
zoom level=6 scale_denom=8735660.37545
zoom level=7 scale_denom=4367830.18772
zoom level=8 scale_denom=2183915.09386
zoom level=9 scale_denom=1091957.54693
zoom level=10 scale_denom=545978.773466
zoom level=11 scale_denom=272989.386733
zoom level=12 scale_denom=136494.693366
zoom level=13 scale_denom=68247.3466832
zoom level=14 scale_denom=34123.6733416
zoom level=15 scale_denom=17061.8366708
zoom level=16 scale_denom=8530.9183354
zoom level=17 scale_denom=4265.4591677
zoom level=18 scale_denom=2132.72958385

I modified tile rendering script to output scales:

#!/opt/python-2_5/bin/python

from math import pi,cos,sin,log,exp,atan

DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi

def minmax (a,b,c):
    a = max(a,b)
    a = min(a,c)
    return a

class GoogleProjection:
    def __init__(self,levels=18):
        self.Bc = []
        self.Cc = []
        self.zc = []
        self.Ac = []
        c = 256
        for d in range(0,levels):
            e = c/2;
            self.Bc.append(c/360.0)
            self.Cc.append(c/(2 * pi))
            self.zc.append((e,e))
            self.Ac.append(c)
            c *= 2
                
    def fromLLtoPixel(self,ll,zoom):
         d = self.zc[zoom]
         e = round(d[0] + ll[0] * self.Bc[zoom])
         f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
         g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
         return (e,g)
     
    def fromPixelToLL(self,px,zoom):
         e = self.zc[zoom]
         f = (px[0] - e[0])/self.Bc[zoom]
         g = (px[1] - e[1])/-self.Cc[zoom]
         h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
         return (f,h)

from mapnik import *

def calculate_scales(minZoom=1,maxZoom=18):
    gprj = GoogleProjection(maxZoom+1)
    proj_merc = Projection("+proj=merc +datum=WGS84")
    
    m = Map(256,256)
    
    for z in range(minZoom,maxZoom + 1):
        px = gprj.fromLLtoPixel((0,0),z)
        x = int(px[0]/256.0)
        y = int(px[1]/256.0)
     
        p0 = gprj.fromPixelToLL((x * 256.0, (y+1) * 256.0),z)
        p1 = gprj.fromPixelToLL(((x+1) * 256.0, y * 256.0),z)
        c0 = proj_merc.forward(Coord(p0[0],p0[1]))
        c1 = proj_merc.forward(Coord(p1[0],p1[1]))
        
        bbox = Envelope(c0.x,c0.y,c1.x,c1.y)
        bbox.width(bbox.width())
        bbox.height(bbox.height())
        m.zoom_to_box(bbox)
        print "zoom level=%s scale_denom=%s " % (z, m.scale()/0.00028)
             
import sys

if __name__ == "__main__":
    argv = sys.argv
    if argv is None or len(argv) != 3:
        print "usage: ./tile_calculator.py minZoom maxZoom "
        sys.exit(0)
        
    minZoom = int(argv[1])
    maxZoom = int(argv[2])
    calculate_scales()
    sys.exit(0)



> Any word 
> on that magical thing you talked about yesterday which would let someone
> interactively play around with styles?l

Well, it is work in progress:). At the moment it can load map definition file 
(e.g osm.xml), navigate around map (zoom in/out/box ,panning), select 
features (see prev link). I'll write about in a sepate thread to see if there 
enough interest to get involved etc etc

Cheers,
--Artem 

>
> * @ 17/01/07 03:23:20 PM steve at asklater.com wrote:
> > * @ 17/01/07 03:21:52 PM steve at asklater.com wrote:
> > > I've added bits so that this layer now auto-renders whatever you look
> > > at providing tiles dont currently exist or are >1 week old. Oh and its
> > > got the latest planet. So look at an area and if its old/not there, its
> > > now in the queue to render. Come back in a bit and it will be. (unless
> > > its a blank tile and its then just thrown away)
> > >
> > > I've switched round the highway colours (I think) so as new tiles
> > > render they get the 'proper' way around.
> >
> > I broke the colours a bit, fixing :)
> >
> > have fun,
> >
> > SteveC steve at asklater.com http://www.asklater.com/steve/
> >
> > _______________________________________________
> > talk mailing list
> > talk at openstreetmap.org
> > http://lists.openstreetmap.org/cgi-bin/mailman/listinfo/talk
>
> have fun,
>
> SteveC steve at asklater.com http://www.asklater.com/steve/