1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

#!/usr/bin/python
#

import mapnik as mapnik
import sys, os

def getShapeFiles(vmdroot,shproot,shpArr = None):
   '''Returns the array of files found in directory tree starting at vmdroot
   whose root file name is shproot.
   '''
   if shpArr == None:
      #print "getShapeFiles - first call : Root"
      shpArr = []

   dirlist = os.listdir(vmdroot)
   for dirent in dirlist:
      dirpath = "%s/%s" % (vmdroot,dirent)
      if os.path.isdir(dirpath):
         #print "%s is a directory - scanning..." % dirpath
         getShapeFiles(dirpath,shproot,shpArr)

      direntsplit = dirent.split('.')
      if direntsplit[0] == shproot:
         if len(direntsplit)>1:
            if direntsplit[1] == 'shp':
               shpArr.append(dirpath)
               sys.stdout.write('.')
               sys.stdout.flush()
               
   return shpArr

def addShpLayer(m,shp,style):
   ''' Add a new layer to mapnik map m using shapefile datasource shp and the
   style named style.
   '''
   shp = shp.split('.')[0]  #strip off file extension
   lyr = mapnik.Layer('shp',"+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +datum=OSGB36 +no_defs no_defs")
   lyr.datasource = mapnik.Shapefile(file='%s' % shp)
   lyr.styles.append(style)
   m.layers.append(lyr)

if __name__ == "__main__":
   mapfile = "styles.xml"
   map_uri = "image.png"

    #---------------------------------------------------
    #  Change this to the bounding box you want
    #
    #ll = (-6.5, 49.5, 2.1, 59)
   Teesside = (-1.3,54.50,-1.17,54.75)
   Hartlepool = (-1.2,54.62,-1.17,54.72)
   ll = Hartlepool
    #---------------------------------------------------

   z = 10
   imgx = 500 * z
   imgy = 500 * z
   
   m = mapnik.Map(imgx,imgy)
   mapnik.load_map(m,mapfile)
   
   #vmdroot = "/home/disk2/os_opendata_vmd/"
   vmdroot = "/home/graham/OSM/data/OS_VMD_NZ"

   print "vmdmap"
   print "Reading Shapefiles"
   # Add Tidal Boundaries
   shpfiles=getShapeFiles(vmdroot,"TidalBoundary")
   for shp in shpfiles:
      addShpLayer(m,shp,'MHW')

    # Add Natural features (water, woodland and other features)
   shpfiles=getShapeFiles(vmdroot,"NaturalFeature_Area")
   for shp in shpfiles:
      addShpLayer(m,shp,'natural-feature-poly')
   shpfiles=getShapeFiles(vmdroot,"NaturalFeature_Line")
   for shp in shpfiles:
      addShpLayer(m,shp,'natural-feature-line')

    # Add Buildings
   shpfiles=getShapeFiles(vmdroot,"Settlement_Area")
   for shp in shpfiles:
      addShpLayer(m,shp,'building')

    # Add Railways
   shpfiles=getShapeFiles(vmdroot,"Railway_Line")
   for shp in shpfiles:
      addShpLayer(m,shp,'railway')
   shpfiles=getShapeFiles(vmdroot,"Railway_Point")
   for shp in shpfiles:
      addShpLayer(m,shp,'railway')

    # Add Roads
   shpfiles=getShapeFiles(vmdroot,"Road_Line")
   for shp in shpfiles:
      addShpLayer(m,shp,'other-road') 
      addShpLayer(m,shp,'secondary')
      addShpLayer(m,shp,'primary')
      addShpLayer(m,shp,'motorway')

    # Add misc text
   shpfiles=getShapeFiles(vmdroot,"Heritage")
   for shp in shpfiles:
      addShpLayer(m,shp,'text')
   shpfiles=getShapeFiles(vmdroot,"Text")
   for shp in shpfiles:
      addShpLayer(m,shp,'text')

   print "."
   print "rendering map..."

   prj = mapnik.Projection("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +datum=OSGB36 +no_defs no_defs")
   c0 = prj.forward(mapnik.Coord(ll[0],ll[1]))
   c1 = prj.forward(mapnik.Coord(ll[2],ll[3]))
   if hasattr(mapnik,'mapnik_version') and mapnik.mapnik_version() >= 800:
      bbox = mapnik.Box2d(c0.x,c0.y,c1.x,c1.y)
   else:
      bbox = mapnik.Envelope(c0.x,c0.y,c1.x,c1.y)
   m.zoom_to_box(bbox)
   im = mapnik.Image(imgx,imgy)
   mapnik.render(m, im)
   view = im.view(0,0,imgx,imgy) # x,y,width,height
   view.save(map_uri,'png')

   print "Done! - output in %s" % map_uri