Форумы GIS-Lab.info

Код: Выделить всё

1	#!/usr/bin/env python
2	from math import pi,cos,sin,log,exp,atan
3	from subprocess import call
4	import sys, os
5	import multiprocessing
6	
7	try:
8	    import mapnik2 as mapnik
9	except:
10	    import mapnik
11	
12	DEG_TO_RAD = pi/180
13	RAD_TO_DEG = 180/pi
14	
15	# Default number of rendering threads to spawn, should be roughly equal to number of CPU cores available
16	NUM_THREADS = 4
17	
18	
19	def minmax (a,b,c):
20	    a = max(a,b)
21	    a = min(a,c)
22	    return a
23	
24	class GoogleProjection:
25	    def __init__(self,levels=18):
26	        self.Bc = []
27	        self.Cc = []
28	        self.zc = []
29	        self.Ac = []
30	        c = 256
31	        for d in range(0,levels):
32	            e = c/2;
33	            self.Bc.append(c/360.0)
34	            self.Cc.append(c/(2 * pi))
35	            self.zc.append((e,e))
36	            self.Ac.append(c)
37	            c *= 2
38	               
39	    def fromLLtoPixel(self,ll,zoom):
40	         d = self.zc[zoom]
41	         e = round(d[0] + ll[0] * self.Bc[zoom])
42	         f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
43	         g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
44	         return (e,g)
45	     
46	    def fromPixelToLL(self,px,zoom):
47	         e = self.zc[zoom]
48	         f = (px[0] - e[0])/self.Bc[zoom]
49	         g = (px[1] - e[1])/-self.Cc[zoom]
50	         h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
51	         return (f,h)
52	
53	
54	
55	class RenderThread:
56	    def __init__(self, tile_dir, mapfile, q, printLock, maxZoom):
57	        self.tile_dir = tile_dir
58	        self.q = q
59	        self.mapfile = mapfile
60	        self.maxZoom = maxZoom
61	        self.printLock = printLock
62	
63	    def render_tile(self, tile_uri, x, y, z):
64	        # Calculate pixel positions of bottom-left & top-right
65	        p0 = (x * 256, (y + 1) * 256)
66	        p1 = ((x + 1) * 256, y * 256)
67	
68	        # Convert to LatLong (EPSG:4326)
69	        l0 = self.tileproj.fromPixelToLL(p0, z);
70	        l1 = self.tileproj.fromPixelToLL(p1, z);
71	
72	        # Convert to map projection (e.g. mercator co-ords EPSG:900913)
73	        c0 = self.prj.forward(mapnik.Coord(l0[0],l0[1]))
74	        c1 = self.prj.forward(mapnik.Coord(l1[0],l1[1]))
75	
76	        # Bounding box for the tile
77	        if hasattr(mapnik,'mapnik_version') and mapnik.mapnik_version() >= 800:
78	            bbox = mapnik.Box2d(c0.x,c0.y, c1.x,c1.y)
79	        else:
80	            bbox = mapnik.Envelope(c0.x,c0.y, c1.x,c1.y)
81	        render_size = 256
82	        self.m.resize(render_size, render_size)
83	        self.m.zoom_to_box(bbox)
84	        self.m.buffer_size = 128
85	
86	        # Render image with default Agg renderer
87	        im = mapnik.Image(render_size, render_size)
88	        mapnik.render(self.m, im)
89	        im.save(tile_uri, 'png256')
90	
91	
92	    def loop(self):
93	       
94	        self.m = mapnik.Map(256, 256)
95	        # Load style XML
96	        mapnik.load_map(self.m, self.mapfile, True)
97	        # Obtain <Map> projection
98	        self.prj = mapnik.Projection(self.m.srs)
99	        # Projects between tile pixel co-ordinates and LatLong (EPSG:4326)
100	        self.tileproj = GoogleProjection(self.maxZoom+1)
101	               
102	        while True:
103	            #Fetch a tile from the queue and render it
104	            r = self.q.get()
105	            if (r == None):
106	                self.q.task_done()
107	                break
108	            else:
109	                (name, tile_uri, x, y, z) = r
110	
111	            exists= ""
112	            if os.path.isfile(tile_uri):
113	                exists= "exists"
114	            else:
115	                self.render_tile(tile_uri, x, y, z)
116	            bytes=os.stat(tile_uri)[6]
117	            empty= ''
118	            if bytes == 103:
119	                empty = " Empty Tile "
120	            self.printLock.acquire()
121	            print name, ":", z, x, y, exists, empty
122	            self.printLock.release()
123	            self.q.task_done()
124	
125	
126	
127	def render_tiles(bbox, mapfile, tile_dir, minZoom=1,maxZoom=18, name="unknown", num_threads=NUM_THREADS):
128	    print "render_tiles(",bbox, mapfile, tile_dir, minZoom,maxZoom, name,")"
129	
130	    # Launch rendering threads
131	    queue = multiprocessing.JoinableQueue(32)
132	    printLock = multiprocessing.Lock()
133	    renderers = {}
134	    for i in range(num_threads):
135	        renderer = RenderThread(tile_dir, mapfile, queue, printLock, maxZoom)
136	        render_thread = multiprocessing.Process(target=renderer.loop)
137	        render_thread.start()
138	        #print "Started render thread %s" % render_thread.getName()
139	        renderers[i] = render_thread
140	
141	    if not os.path.isdir(tile_dir):
142	         os.mkdir(tile_dir)
143	
144	    gprj = GoogleProjection(maxZoom+1) 
145	
146	    ll0 = (bbox[0],bbox[3])
147	    ll1 = (bbox[2],bbox[1])
148	
149	    for z in range(minZoom,maxZoom + 1):
150	        px0 = gprj.fromLLtoPixel(ll0,z)
151	        px1 = gprj.fromLLtoPixel(ll1,z)
152	
153	        # check if we have directories in place
154	        zoom = "%s" % z
155	        if not os.path.isdir(tile_dir + zoom):
156	            os.mkdir(tile_dir + zoom)
157	        for x in range(int(px0[0]/256.0),int(px1[0]/256.0)+1):
158	            # Validate x co-ordinate
159	            if (x < 0) or (x >= 2**z):
160	                continue
161	            # check if we have directories in place
162	            str_x = "%s" % x
163	            if not os.path.isdir(tile_dir + zoom + '/' + str_x):
164	                os.mkdir(tile_dir + zoom + '/' + str_x)
165	            for y in range(int(px0[1]/256.0),int(px1[1]/256.0)+1):
166	                # Validate x co-ordinate
167	                if (y < 0) or (y >= 2**z):
168	                    continue
169	                str_y = "%s" % y
170	                tile_uri = tile_dir + zoom + '/' + str_x + '/' + str_y + '.png'
171	                # Submit tile to be rendered into the queue
172	                t = (name, tile_uri, x, y, z)
173	                queue.put(t)
174	
175	    # Signal render threads to exit by sending empty request to queue
176	    for i in range(num_threads):
177	        queue.put(None)
178	    # wait for pending rendering jobs to complete
179	    queue.join()
180	    for i in range(num_threads):
181	        renderers[i].join()
182	
183	
184	
185	if __name__ == "__main__":
186	       
187	    home = os.environ['HOME']
188	    try:
189	        mapfile = os.environ['MAPNIK_MAP_FILE']
190	    except KeyError:
191	        mapfile = home + "/svn.openstreetmap.org/applications/rendering/mapnik/osm-local.xml"
192	    try:
193	        tile_dir = os.environ['MAPNIK_TILE_DIR']
194	    except KeyError:
195	        tile_dir = home + "/osm/tiles/"
196	
197	    if not tile_dir.endswith('/'):
198	        tile_dir = tile_dir + '/'
199	
200	    #-------------------------------------------------------------------------
201	    #
202	    # Change the following for different bounding boxes and zoom levels
203	    #
204	    # Start with an overview
205	    # World
206	    bbox = (-180.0,-90.0, 180.0,90.0)
207	
208	    render_tiles(bbox, mapfile, tile_dir, 0, 5, "World")
209	
210	    minZoom = 10
211	    maxZoom = 16
212	    bbox = (-2, 50.0,1.0,52.0)
213	    render_tiles(bbox, mapfile, tile_dir, minZoom, maxZoom)
214	
215	    # Muenchen
216	    bbox = (11.4,48.07, 11.7,48.22)
217	    render_tiles(bbox, mapfile, tile_dir, 1, 12 , "Muenchen")
218	
219	    # Muenchen+
220	    bbox = (11.3,48.01, 12.15,48.44)
221	    render_tiles(bbox, mapfile, tile_dir, 7, 12 , "Muenchen+")
222	
223	    # Muenchen++
224	    bbox = (10.92,47.7, 12.24,48.61)
225	    render_tiles(bbox, mapfile, tile_dir, 7, 12 , "Muenchen++")
226	
227	    # Nuernberg
228	    bbox=(10.903198,49.560441,49.633534,11.038085)
229	    render_tiles(bbox, mapfile, tile_dir, 10, 16, "Nuernberg")
230	
231	    # Karlsruhe
232	    bbox=(8.179113,48.933617,8.489252,49.081707)
233	    render_tiles(bbox, mapfile, tile_dir, 10, 16, "Karlsruhe")
234	
235	    # Karlsruhe+
236	    bbox = (8.3,48.95,8.5,49.05)
237	    render_tiles(bbox, mapfile, tile_dir, 1, 16, "Karlsruhe+")
238	
239	    # Augsburg
240	    bbox = (8.3,48.95,8.5,49.05)
241	    render_tiles(bbox, mapfile, tile_dir, 1, 16, "Augsburg")
242	
243	    # Augsburg+
244	    bbox=(10.773251,48.369594,10.883834,48.438577)
245	    render_tiles(bbox, mapfile, tile_dir, 10, 14, "Augsburg+")
246	
247	    # Europe+
248	    bbox = (1.0,10.0, 20.6,50.0)
249	    render_tiles(bbox, mapfile, tile_dir, 1, 11 , "Europe+")