1 /***************************************
2 Data extractor and file dumper.
3
4 Part of the Routino routing software.
5 ******************/ /******************
6 This file Copyright 2008-2017 Andrew M. Bishop
7
8 This program is free software: you can redistribute it and/or modify
9 it under the terms of the GNU Affero General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU Affero General Public License for more details.
17
18 You should have received a copy of the GNU Affero General Public License
19 along with this program. If not, see <http://www.gnu.org/licenses/>.
20 ***************************************/
21
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include <sys/stat.h>
27 #include <sys/time.h>
28 #include <time.h>
29 #include <math.h>
30
31 #include "types.h"
32 #include "functions.h"
33 #include "nodes.h"
34 #include "segments.h"
35 #include "ways.h"
36 #include "sorting.h"
37
38
39 /* Global variables (required to link sorting.c) */
40
41 /*+ The command line '--tmpdir' option or its default value. +*/
42 char *option_tmpdirname=NULL;
43
44 /*+ The amount of RAM to use for filesorting. +*/
45 size_t option_filesort_ramsize=0;
46
47 /*+ The number of threads to use for filesorting. +*/
48 int option_filesort_threads=1;
49
50 /* Local types */
51
52 typedef struct _crossing
53 {
54 double f;
55 uint32_t x;
56 uint32_t y;
57 }
58 crossing;
59
60 /* Local functions */
61
62 static void ll_xy_rad(int zoom,double longitude,double latitude,double *x,double *y);
63
64 static int compare_crossing(crossing *a,crossing *b);
65
66
67 /*++++++++++++++++++++++++++++++++++++++
68 The main function.
69 ++++++++++++++++++++++++++++++++++++++*/
70
main(int argc,char ** argv)71 int main(int argc,char** argv)
72 {
73 Nodes *OSMNodes;
74 Segments *OSMSegments;
75 Ways *OSMWays;
76 int arg;
77 char *dirname=NULL,*prefix=NULL;
78 char *nodes_filename,*segments_filename,*ways_filename;
79 index_t item;
80 crossing *crossings=(crossing*)malloc(128*sizeof(crossing));
81 crossing **crossingsp=(crossing**)malloc(128*sizeof(crossing*));
82 int crossing_alloc=128;
83 float ***highways,***transports,***properties,***speeds;
84 double lat,lon,x,y;
85 uint32_t xmin,ymin,xmax,ymax,xt,yt;
86 int zoom=13;
87 int speed_count=0;
88 double *speed_selection=NULL;
89
90 /* Parse the command line arguments */
91
92 for(arg=1;arg<argc;arg++)
93 {
94 if(!strcmp(argv[arg],"--help"))
95 goto usage;
96 else if(!strncmp(argv[arg],"--dir=",6))
97 dirname=&argv[arg][6];
98 else if(!strncmp(argv[arg],"--prefix=",9))
99 prefix=&argv[arg][9];
100 else if(!strncmp(argv[arg],"--speed=",8))
101 {
102 int i;
103
104 speed_count=1;
105
106 for(i=8;argv[arg][i];i++)
107 if(argv[arg][i]==',')
108 speed_count++;
109
110 speed_selection=calloc(speed_count,sizeof(double));
111
112 speed_count=0;
113
114 for(i=7;argv[arg][i];i++)
115 if(argv[arg][i]==',' || argv[arg][i]=='=')
116 speed_selection[speed_count++]=atof(&argv[arg][i+1]);
117 }
118 else if(!strncmp(argv[arg],"--zoom=",7))
119 {
120 zoom=atoi(&argv[arg][7]);
121
122 if(zoom<10 || zoom>16)
123 goto usage;
124 }
125 else
126 {
127 usage:
128
129 fprintf(stderr,"Usage: dumper\n"
130 " [--help]\n"
131 " [--dir=<name>] [--prefix=<name>]\n"
132 " [--speed=<speed1>,<speed2>,...,<speedN>]"
133 " [--zoom=<10-16>]\n");
134
135 return(1);
136 }
137 }
138
139 /* Load in the data */
140
141 OSMNodes=LoadNodeList(nodes_filename=FileName(dirname,prefix,"nodes.mem"));
142
143 if(!OSMNodes)
144 {
145 fprintf(stderr,"Cannot open nodes file '%s'.\n",nodes_filename);
146 return(1);
147 }
148
149 OSMSegments=LoadSegmentList(segments_filename=FileName(dirname,prefix,"segments.mem"));
150
151 if(!OSMSegments)
152 {
153 fprintf(stderr,"Cannot open segments file '%s'.\n",segments_filename);
154 return(1);
155 }
156
157 OSMWays=LoadWayList(ways_filename=FileName(dirname,prefix,"ways.mem"));
158
159 if(!OSMWays)
160 {
161 fprintf(stderr,"Cannot open ways file '%s'.\n",ways_filename);
162 return(1);
163 }
164
165 /* Allocate the arrays */
166
167 lat=latlong_to_radians(bin_to_latlong(OSMNodes->file.latzero));
168 lon=latlong_to_radians(bin_to_latlong(OSMNodes->file.lonzero));
169
170 ll_xy_rad(zoom,lon,lat,&x,&y);
171
172 xmin=(uint32_t)floor(x);
173 ymax=(uint32_t)floor(y);
174
175 lat=latlong_to_radians(bin_to_latlong(OSMNodes->file.latzero+OSMNodes->file.latbins));
176 lon=latlong_to_radians(bin_to_latlong(OSMNodes->file.lonzero+OSMNodes->file.lonbins));
177
178 ll_xy_rad(zoom,lon,lat,&x,&y);
179
180 xmax=(uint32_t)floor(x);
181 ymin=(uint32_t)floor(y);
182
183 highways =malloc(sizeof(highways[0]) *(xmax-xmin+1));
184 transports=malloc(sizeof(transports[0])*(xmax-xmin+1));
185 properties=malloc(sizeof(properties[0])*(xmax-xmin+1));
186 speeds =malloc(sizeof(speeds[0]) *(xmax-xmin+1));
187
188 for(xt=xmin;xt<=xmax;xt++)
189 {
190 highways[xt-xmin] =malloc(sizeof(highways[0][0]) *(ymax-ymin+1));
191 transports[xt-xmin]=malloc(sizeof(transports[0][0])*(ymax-ymin+1));
192 properties[xt-xmin]=malloc(sizeof(properties[0][0])*(ymax-ymin+1));
193 speeds[xt-xmin] =malloc(sizeof(speeds[0][0]) *(ymax-ymin+1));
194
195 for(yt=ymin;yt<=ymax;yt++)
196 {
197 highways[xt-xmin][yt-ymin] =calloc((Highway_Count-1) ,sizeof(highways[0][0][0]) );
198 transports[xt-xmin][yt-ymin]=calloc((Transport_Count-1),sizeof(transports[0][0][0]));
199 properties[xt-xmin][yt-ymin]=calloc((Property_Count+1) ,sizeof(properties[0][0][0]));
200 speeds[xt-xmin][yt-ymin] =calloc((speed_count) ,sizeof(speeds[0][0][0]) );
201 }
202 }
203
204 /* Dump the segments out with lat/long in tile units. */
205
206 for(item=0;item<OSMSegments->file.number;item++)
207 {
208 Segment *segment=LookupSegment(OSMSegments,item,1);
209
210 if(IsNormalSegment(segment))
211 {
212 double latitude1,longitude1,latitude2,longitude2;
213 double x1,y1,x2,y2;
214 uint32_t xi1,yi1,xi2,yi2;
215 Way *way;
216 double distance;
217
218 /* Get the coordinates */
219
220 distance=1000*distance_to_km(DISTANCE(segment->distance));
221
222 GetLatLong(OSMNodes,segment->node1,NULL,&latitude1,&longitude1);
223
224 GetLatLong(OSMNodes,segment->node2,NULL,&latitude2,&longitude2);
225
226 way=LookupWay(OSMWays,segment->way,1);
227
228 ll_xy_rad(zoom,longitude1,latitude1,&x1,&y1);
229 ll_xy_rad(zoom,longitude2,latitude2,&x2,&y2);
230
231 /* Map to tiles and find tile crossings */
232
233 xi1=(uint32_t)floor(x1);
234 yi1=(uint32_t)floor(y1);
235
236 xi2=(uint32_t)floor(x2);
237 yi2=(uint32_t)floor(y2);
238
239 if(xi1==xi2 && yi1==yi2)
240 {
241 int j;
242
243 for(j=1;j<Highway_Count;j++)
244 if(HIGHWAY(way->type) == j)
245 highways[xi1-xmin][yi1-ymin][j-1]+=distance;
246
247 for(j=1;j<Transport_Count;j++)
248 if(way->allow & TRANSPORTS(j))
249 transports[xi1-xmin][yi1-ymin][j-1]+=distance;
250
251 for(j=1;j<Property_Count;j++)
252 if(way->props & PROPERTIES(j))
253 properties[xi1-xmin][yi1-ymin][j-1]+=distance;
254
255 if(way->type & Highway_CycleBothWays)
256 properties[xi1-xmin][yi1-ymin][Property_Count-1]+=distance;
257
258 if(way->type & Highway_OneWay)
259 properties[xi1-xmin][yi1-ymin][Property_Count]+=distance;
260
261 for(j=0;j<speed_count;j++)
262 if(speed_to_kph(way->speed) == speed_selection[j])
263 speeds[xi1-xmin][yi1-ymin][j]+=distance;
264 }
265 else
266 {
267 int crossing_count=2,i;
268 double lastf;
269 uint32_t lastx,lasty;
270
271 crossings[0].f=0.0;
272 crossings[0].x=xi1;
273 crossings[0].y=yi1;
274
275 crossings[1].f=1.0;
276 crossings[1].x=xi2;
277 crossings[1].y=yi2;
278
279 /* Find X boundaries */
280
281 if(xi1!=xi2)
282 {
283 uint32_t x,minx,maxx;
284
285 if(xi1<xi2)
286 {
287 minx=xi1+1;
288 maxx=xi2;
289 }
290 else
291 {
292 minx=xi2+1;
293 maxx=xi1;
294 }
295
296 for(x=minx;x<=maxx;x++)
297 {
298 double fraction=(double)(x-x1)/(double)(x2-x1);
299 uint32_t y=(uint32_t)(y1+fraction*(y2-y1));
300
301 crossings[crossing_count].f=fraction;
302 crossings[crossing_count].x=x;
303 crossings[crossing_count].y=y;
304 crossing_count++;
305
306 if(crossing_count==crossing_alloc)
307 {
308 crossing_alloc+=128;
309 crossings =(crossing *)realloc(crossings ,crossing_alloc*sizeof(crossing));
310 crossingsp=(crossing**)realloc(crossingsp,crossing_alloc*sizeof(crossing*));
311 }
312 }
313 }
314
315 /* Find Y boundaries */
316
317 if(yi1!=yi2)
318 {
319 uint32_t y,miny,maxy;
320
321 if(yi1<yi2)
322 {
323 miny=yi1+1;
324 maxy=yi2;
325 }
326 else
327 {
328 miny=yi2+1;
329 maxy=yi1;
330 }
331
332 for(y=miny;y<=maxy;y++)
333 {
334 double fraction=(double)(y-y1)/(double)(y2-y1);
335 uint32_t x=(uint32_t)(x1+fraction*(x2-x1));
336
337 crossings[crossing_count].f=fraction;
338 crossings[crossing_count].x=x;
339 crossings[crossing_count].y=y;
340 crossing_count++;
341
342 if(crossing_count==crossing_alloc)
343 {
344 crossing_alloc+=128;
345 crossings =(crossing *)realloc(crossings ,crossing_alloc*sizeof(crossing));
346 crossingsp=(crossing**)realloc(crossingsp,crossing_alloc*sizeof(crossing*));
347 }
348 }
349 }
350
351 /* Sort the boundaries */
352
353 for(i=0;i<crossing_count;i++)
354 crossingsp[i]=&crossings[i];
355
356 filesort_heapsort((void**)crossingsp,crossing_count,(int (*)(const void*,const void*))compare_crossing);
357
358 /* Partition the data among the tiles */
359
360 lastf=crossingsp[0]->f;
361 lastx=crossingsp[0]->x;
362 lasty=crossingsp[0]->y;
363
364 for(i=1;i<crossing_count;i++)
365 {
366 double f=crossingsp[i]->f;
367 uint32_t xi,x=crossingsp[i]->x;
368 uint32_t yi,y=crossingsp[i]->y;
369 double d;
370 int j;
371
372 xi=x;
373 if(lastx<x) xi=lastx;
374
375 yi=y;
376 if(lasty<y) yi=lasty;
377
378 d=distance*(f-lastf);
379
380 for(j=1;j<Highway_Count;j++)
381 if(HIGHWAY(way->type) == j)
382 highways[xi-xmin][yi-ymin][j-1]+=d;
383
384 for(j=1;j<Transport_Count;j++)
385 if(way->allow & TRANSPORTS(j))
386 transports[xi-xmin][yi-ymin][j-1]+=d;
387
388 for(j=1;j<Property_Count;j++)
389 if(way->props & PROPERTIES(j))
390 properties[xi-xmin][yi-ymin][j-1]+=d;
391
392 if(way->type & Highway_CycleBothWays)
393 properties[xi1-xmin][yi1-ymin][Property_Count-1]+=distance;
394
395 if(way->type & Highway_OneWay)
396 properties[xi1-xmin][yi1-ymin][Property_Count]+=distance;
397
398 for(j=0;j<speed_count;j++)
399 if(speed_to_kph(way->speed) == speed_selection[j])
400 speeds[xi1-xmin][yi1-ymin][j]+=distance;
401
402 lastf=f;
403 lastx=x;
404 lasty=y;
405 }
406 }
407 }
408 }
409
410 /* Print the results */
411
412 for(xt=xmin;xt<=xmax;xt++)
413 {
414 for(yt=ymin;yt<=ymax;yt++)
415 {
416 int do_highways=0,do_transports=0,do_properties=0;
417 int j;
418
419 for(j=1;j<Highway_Count;j++)
420 if(highways[xt-xmin][yt-ymin][j-1] > 0)
421 do_highways++;
422
423 for(j=1;j<Transport_Count;j++)
424 if(transports[xt-xmin][yt-ymin][j-1] > 0)
425 do_transports++;
426
427 for(j=1;j<Property_Count;j++)
428 if(properties[xt-xmin][yt-ymin][j-1] > 0)
429 do_properties++;
430
431 if(do_highways || do_transports || do_properties)
432 {
433 printf("%d %d",xt,yt);
434
435 for(j=1;j<Highway_Count;j++)
436 printf(" %.2f",highways[xt-xmin][yt-ymin][j-1]);
437
438 for(j=1;j<Transport_Count;j++)
439 printf(" %.2f",transports[xt-xmin][yt-ymin][j-1]);
440
441 for(j=1;j<=Property_Count+1;j++)
442 printf(" %.2f",properties[xt-xmin][yt-ymin][j-1]);
443
444 for(j=0;j<speed_count;j++)
445 printf(" %.2f",speeds[xt-xmin][yt-ymin][j]);
446
447 printf("\n");
448 }
449
450 free(highways[xt-xmin][yt-ymin]);
451 free(transports[xt-xmin][yt-ymin]);
452 free(properties[xt-xmin][yt-ymin]);
453 free(speeds[xt-xmin][yt-ymin]);
454 }
455
456 free(highways[xt-xmin]);
457 free(transports[xt-xmin]);
458 free(properties[xt-xmin]);
459 free(speeds[xt-xmin]);
460 }
461
462 free(highways);
463 free(transports);
464 free(properties);
465 free(speeds);
466
467 free(speed_selection);
468
469 free(crossings);
470 free(crossingsp);
471
472 return(0);
473 }
474
475
476 /*++++++++++++++++++++++++++++++++++++++
477 Sort the tile edge crossings into order.
478
479 int compare_crossing Returns the comparison of the fraction (f) field.
480
481 crossing *a The first boundary crossing.
482
483 crossing *b The second boundary crossing.
484 ++++++++++++++++++++++++++++++++++++++*/
485
compare_crossing(crossing * a,crossing * b)486 static int compare_crossing(crossing *a,crossing *b)
487 {
488 double ad=a->f;
489 double bd=b->f;
490
491 if(ad<bd)
492 return(-1);
493 else if(ad>bd)
494 return(1);
495 else
496 return(-FILESORT_PRESERVE_ORDER(a,b)); /* latest version first */
497 }
498
499
500 /*++++++++++++++++++++++++++++++++++++++
501 Convert latitude and longitude into tile coordinates.
502
503 int zoom The zoom level.
504
505 double longitude The latitude of the point.
506
507 double latitude The longitude of the point.
508
509 double *x The tile X number (including fractional part).
510
511 double *y The tile Y number (including fractional part).
512 ++++++++++++++++++++++++++++++++++++++*/
513
ll_xy_rad(int zoom,double longitude,double latitude,double * x,double * y)514 static void ll_xy_rad(int zoom,double longitude,double latitude,double *x,double *y)
515 {
516 *x=longitude/M_PI+1;
517 *y=1-log(tan(latitude/2+(M_PI/4)))/M_PI;
518
519 *x/=pow(2,1-zoom);
520 *y/=pow(2,1-zoom);
521 }
522