1 | # -*- coding: utf-8 -*- |
---|
2 | # by kay - basic functions |
---|
3 | |
---|
4 | ### config for Trunk |
---|
5 | DSN = 'dbname=gis' |
---|
6 | |
---|
7 | import sys |
---|
8 | import commands |
---|
9 | import psycopg2 |
---|
10 | import csv |
---|
11 | import re |
---|
12 | from numpy import * |
---|
13 | |
---|
14 | highwaytypes = { |
---|
15 | 'trunk':['<0.9,1,0.9>',1.0], |
---|
16 | 'trunk_link':['<0.9,1,0.9>',1.0], |
---|
17 | 'primary':['<1,1,0.9>',1.0], |
---|
18 | 'primary_link':['<1,1,0.9>',1.0], |
---|
19 | 'secondary':['<1,0.9,0.9>',1.0], |
---|
20 | 'secondary_link':['<1,0.9,0.9>',1.0], |
---|
21 | 'tertiary':['<1,0.9,0.8>',1.0], |
---|
22 | 'residential':['<0.9,0.9,0.9>',0.8], |
---|
23 | 'living_street':['<0.8,0.8,0.9>',0.8], |
---|
24 | 'unclassified':['<0.8,0.8,0.8>',0.8], |
---|
25 | 'service':['<0.8,0.8,0.8>',0.6], |
---|
26 | 'pedestrian':['<0.8,0.9,0.8>',0.8], |
---|
27 | 'footway':['<0.8,0.9,0.8>',0.3], |
---|
28 | 'cycleway':['<0.8,0.8,0.9>',0.3], |
---|
29 | 'path':['<0.8,0.9,0.8>',0.5] |
---|
30 | } |
---|
31 | buildingtypes = { |
---|
32 | 'yes':['<1,0.8,0.8>',1.0], |
---|
33 | 'office':['<1,0.6,0.6>',1.0], |
---|
34 | 'apartments':['<0.8,1,0.6>',1.0], |
---|
35 | 'garages':['<0.8,0.8,0.8>',1.0], |
---|
36 | } |
---|
37 | amenitybuildingtypes = { |
---|
38 | 'place_of_worship':['<1,1,0.6>',1.0], |
---|
39 | 'hospital':['<1,0.6,0.6>',1.0], |
---|
40 | 'theatre':['<1,0.6,1>',1.0], |
---|
41 | 'university':['<0.6,1,0.8>',1.0], |
---|
42 | 'parking':['<0.6,0.6,1>',1.0] |
---|
43 | } |
---|
44 | |
---|
45 | def avg(a,b): return (a+b)/2.0 |
---|
46 | |
---|
47 | def shift_by_meters(lat, lon, brng, d): |
---|
48 | R = 6371000.0 # earth's radius in m |
---|
49 | lat=math.radians(lat) |
---|
50 | lon=math.radians(lon) |
---|
51 | lat2 = math.asin( math.sin(lat)*math.cos(d/R) + |
---|
52 | math.cos(lat)*math.sin(d/R)*math.cos(brng)) |
---|
53 | lon2 = lon + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat), |
---|
54 | math.cos(d/R)-math.sin(lat)*math.sin(lat2)) |
---|
55 | lat2=math.degrees(lat2) |
---|
56 | lon2=math.degrees(lon2) |
---|
57 | return [lat2,lon2] |
---|
58 | |
---|
59 | def calc_bearing(x1,y1,x2,y2,side): |
---|
60 | Q = complex(x1,y1) |
---|
61 | R = complex(x2,y2) |
---|
62 | v = R-Q |
---|
63 | if side=='left': |
---|
64 | v=v*complex(0,1); |
---|
65 | elif side=='right': |
---|
66 | v=v*complex(0,-1); |
---|
67 | else: |
---|
68 | raise TypeError('side must be left or right') |
---|
69 | #v=v*(1/abs(v)) # normalize |
---|
70 | angl = angle(v) # angle (radians) (0°=right, and counterclockwise) |
---|
71 | bearing = math.pi/2.0-angl # (0°=up, and clockwise) |
---|
72 | return bearing |
---|
73 | |
---|
74 | def pov_highway(f,highway): |
---|
75 | highwaytype = highway[1] |
---|
76 | #highwaytype = 'secondary' |
---|
77 | highwayparams = highwaytypes.get(highwaytype) |
---|
78 | linestring = highway[2] |
---|
79 | linestring = linestring[11:] # cut off the "LINESTRING(" |
---|
80 | linestring = linestring[:-1] # cut off the ")" |
---|
81 | points = linestring.split(',') |
---|
82 | |
---|
83 | oneway = False |
---|
84 | if(highway[5]=='yes'): |
---|
85 | oneway=True |
---|
86 | |
---|
87 | lanes = highway[3] |
---|
88 | if lanes==None: |
---|
89 | lanefactor=2.0 # 2 lanes seems to be default |
---|
90 | if oneway: |
---|
91 | lanefactor=1.2 # except for one-ways (make them a bit wider than 2 lanes/2) |
---|
92 | else: |
---|
93 | lanefactor=float(lanes) |
---|
94 | lanewidth = 2.5 # m |
---|
95 | streetwidth = lanewidth * lanefactor |
---|
96 | |
---|
97 | layer = highway[4] |
---|
98 | if layer==None: |
---|
99 | layer='0' |
---|
100 | layer = int(layer) |
---|
101 | if layer<0: |
---|
102 | layer=0 # FIXME |
---|
103 | layerheight = 4.0*layer # 4 m per layer |
---|
104 | layerheight = layerheight / 0.05 # counteract the scale statement |
---|
105 | |
---|
106 | numpoints = len(points) |
---|
107 | |
---|
108 | # draw road |
---|
109 | f.write("sphere_sweep {{ linear_spline, {0},\n".format(numpoints+2)) |
---|
110 | f.write("/* osm_id={0} */\n".format(highway[0])) |
---|
111 | |
---|
112 | for i,point in enumerate(points): |
---|
113 | latlon = point.split(' ') |
---|
114 | if (i==0): |
---|
115 | f.write(" <{0}, {3}, {1}>,{2}\n".format(latlon[0],latlon[1],streetwidth*highwayparams[1],layerheight)) |
---|
116 | f.write(" <{0}, {3}, {1}>,{2}\n".format(latlon[0],latlon[1],streetwidth*highwayparams[1],layerheight)) |
---|
117 | if (i==numpoints-1): |
---|
118 | f.write(" <{0}, {3}, {1}>,{2}\n".format(latlon[0],latlon[1],streetwidth*highwayparams[1],layerheight)) |
---|
119 | |
---|
120 | print highwayparams[0],highwayparams[1],streetwidth |
---|
121 | f.write(""" tolerance 1 |
---|
122 | texture {{ |
---|
123 | pigment {{ |
---|
124 | color rgb {0} |
---|
125 | }} |
---|
126 | finish {{ |
---|
127 | specular 0.05 |
---|
128 | roughness 0.05 |
---|
129 | /*reflection 0.5*/ |
---|
130 | }} |
---|
131 | }} |
---|
132 | scale <1, 0.05, 1> |
---|
133 | }} |
---|
134 | \n""".format(highwayparams[0])) |
---|
135 | # draw casing |
---|
136 | f.write("sphere_sweep {{ linear_spline, {0},\n".format(numpoints+2)) |
---|
137 | f.write("/* osm_id={0} */\n".format(highway[0])) |
---|
138 | |
---|
139 | for i,point in enumerate(points): |
---|
140 | latlon = point.split(' ') |
---|
141 | if (i==0): |
---|
142 | f.write(" <{0}, {3}, {1}>,{2}\n".format(latlon[0],latlon[1],1.2*streetwidth*highwayparams[1],layerheight)) |
---|
143 | f.write(" <{0}, {3}, {1}>,{2}\n".format(latlon[0],latlon[1],1.2*streetwidth*highwayparams[1],layerheight)) |
---|
144 | if (i==numpoints-1): |
---|
145 | f.write(" <{0}, {3}, {1}>,{2}\n".format(latlon[0],latlon[1],1.2*streetwidth*highwayparams[1],layerheight)) |
---|
146 | |
---|
147 | print highwayparams[0],highwayparams[1],streetwidth |
---|
148 | f.write(""" tolerance 1 |
---|
149 | texture {{ |
---|
150 | pigment {{ |
---|
151 | color rgb <0.2,0.2,0.2> |
---|
152 | }} |
---|
153 | finish {{ |
---|
154 | specular 0.05 |
---|
155 | roughness 0.05 |
---|
156 | /*reflection 0.5*/ |
---|
157 | }} |
---|
158 | }} |
---|
159 | scale <1, 0.05, 1> |
---|
160 | translate <0, -0.1, 0> |
---|
161 | }} |
---|
162 | \n""".format(highwayparams[0])) |
---|
163 | |
---|
164 | def parse_length_in_meters(length,default): |
---|
165 | units={ |
---|
166 | 'm':1.0, 'metres':1.0, 'meters':1.0, 'metre':1.0, 'meter':1.0, |
---|
167 | 'km':1000.0, |
---|
168 | 'ft':0.3, |
---|
169 | 'yd':1.1, |
---|
170 | } |
---|
171 | parsed = re.split('(\d*[,.]?\d+)',length) |
---|
172 | if len(parsed)!=3: |
---|
173 | print "### unparsable length '{0}', parsed: {1}".format(length,parsed) |
---|
174 | return default |
---|
175 | prefix = parsed[0].strip() |
---|
176 | if prefix!='': |
---|
177 | print "### unknown prefix '{0}' in length '{1}'".format(prefix,length) |
---|
178 | return default |
---|
179 | unit = parsed[2].strip().lower() |
---|
180 | if unit=='': unit='m' # defaults to m |
---|
181 | factor = units.get(unit,0.0) |
---|
182 | if factor==0.0: |
---|
183 | print "### unknown unit '{0}' in length '{1}'".format(unit,length) |
---|
184 | meter = float(parsed[1])*factor |
---|
185 | return meter |
---|
186 | |
---|
187 | def pov_building(f,building): |
---|
188 | buildingtype = building[1] |
---|
189 | #buildingtype = 'secondary' |
---|
190 | buildingparams = buildingtypes.get(buildingtype) |
---|
191 | |
---|
192 | linestring = building[2] |
---|
193 | linestring = linestring[8:] # cut off the "POLYGON(" |
---|
194 | linestring = linestring[:-1] # cut off the ")" |
---|
195 | |
---|
196 | heightstring = building[3] |
---|
197 | print 'hs="{0}"'.format(heightstring) |
---|
198 | if (heightstring==None): |
---|
199 | height = 10.0 |
---|
200 | else: |
---|
201 | height = parse_length_in_meters(heightstring,0.01) |
---|
202 | |
---|
203 | amenity = building[4] |
---|
204 | print amenity |
---|
205 | |
---|
206 | points = linestring.split(',')#.strip('(').strip(')') |
---|
207 | #print points |
---|
208 | |
---|
209 | numpoints = len(points) |
---|
210 | f.write("prism {{ linear_spline 0, 1, {0},\n".format(numpoints)) |
---|
211 | f.write("/* osm_id={0} */\n".format(building[0])) |
---|
212 | |
---|
213 | for i,point in enumerate(points): |
---|
214 | latlon = point.split(' ') |
---|
215 | if (i==0): |
---|
216 | firstpoint="<{0}, {1}>\n".format(latlon[0],latlon[1]) |
---|
217 | if (i!=numpoints-1): |
---|
218 | f.write(" <{0}, {1}>,\n".format(latlon[0].strip('(').strip(')'),latlon[1].strip('(').strip(')'))) |
---|
219 | else: |
---|
220 | f.write(" <{0}, {1}>\n".format(latlon[0].strip('(').strip(')'),latlon[1].strip('(').strip(')'))) |
---|
221 | #f.write(firstpoint) |
---|
222 | |
---|
223 | color = buildingparams[0] |
---|
224 | if amenitybuildingtypes.has_key(amenity): # if amenity is known, use that color |
---|
225 | amenitybuildingparams = amenitybuildingtypes.get(amenity) |
---|
226 | color = amenitybuildingparams[0] |
---|
227 | |
---|
228 | if height != 10.0: |
---|
229 | print 'height:', height |
---|
230 | #color = '<0,1,0>' |
---|
231 | f.write(""" |
---|
232 | texture {{ |
---|
233 | pigment {{ |
---|
234 | color rgb {0} |
---|
235 | }} |
---|
236 | finish {{ |
---|
237 | specular 0.5 |
---|
238 | roughness 0.05 |
---|
239 | /*reflection 0.5*/ |
---|
240 | }} |
---|
241 | }} |
---|
242 | scale <1, {1}, 1> |
---|
243 | }} |
---|
244 | \n""".format(color,height)) |
---|
245 | |
---|