git.maemo.org Git - navit-package/blob - navit/script/osm2navit_sea/osm2navit_sea.pl

   1 #!usr/bin/perl -w
   2 #
   3 # Vincent "latouche" Touchard, vtouchar@latouche.info
   4 # 2008/09/10
   5 #
   6 # Transform a .osm file adding polygons representing the sea and islands.
   7 # The sea is a way with the tag natural=coastline_outer and the islands are
   8 # ways with the tag natural=coastline_inner.
   9 # Adding those tags to osm2navit allow the see to be viewed by navit
  10 #
  11 # To avoid big polygons, the map is divided into small tiles (0.5°x0.5°)
  12 #
  13 # usage:
  14 #   osm2navit_sea.pl minlat minlon maxlat maxlon source_file output_file
  15 #
  16 # lisense: GPL
  17 #
  18 # history:
  19 #  v1: first release
  20 #  v2: 2008/09/06
  21 #      - map split into tiles of 0.5x0.5°
  22 #  v3: 2008/09/10
  23 #      - detects islands vs lakes
  24 #        takes each polygon and try to find is it is inside an other one. If true,
  25 #        it is an island, if not, it is the sea. Should be a little more clever to detect
  26 #        lake (tagged as natural=coastline instead of natural=water) inside islands.
  27 #        This algo is not perfect but works in most cases
  28 #      - sea tile with no way crossing are now drawn in a more clever way:
  29 #        a node for each corner is added to nodes arrays and a way linking those
  30 #        node is added. Those components are then added to the output file during
  31 #        the last section of this script. This allow to detects islands in such tiles
  32 #        as islands and not sea
  33 #  v4: 2009/03/14
  34 #      - new way to compute the intersection of a segment and the border of a tile
  35 #      - don't use nodes that are out of the bbox specified as argument
  36 #      - use lazy regex to parse the .osm file
  37 #
  38 # TODO:
  39 #  - better lookup_handler (cf this fuction for more info)
  40 #  - detect clockwise/anticlockwise ways => sea or land
  41 #     for now -> every closed area is considered as land which is false.
  42 #     some lakes are tags as natural=coastline instead of natural=water
  43 #     ==> done, the algo should be improved for better results
  44
  45 use strict;
  46 use warnings;
  47
  48 use Math::Trig;
  49 use POSIX qw(ceil);
  50 use way;
  51
  52 # lon € [-180, 180]
  53 # lat € [-85.0511, 85.0511]
  54 if ( @ARGV ne 6 ||
  55         abs($ARGV[0]) > 85.0511 || abs($ARGV[2]) > 85.0511 ||
  56         abs($ARGV[1]) > 180 || abs($ARGV[3]) > 180) {
  57         print "Usage: osm2navit_sea.pl minlat minlon maxlat maxlon source_file output_file\n";
  58         exit(1);
  59 }
  60
  61
  62 my $source=$ARGV[4];
  63 my $destination=$ARGV[5];
  64
  65 my ($minlat, $minlong, $maxlat, $maxlong);
  66 my ($minlat2, $minlong2, $maxlat2, $maxlong2);
  67
  68 $minlat=$ARGV[0];
  69 $minlong=$ARGV[1];
  70 $maxlat=$ARGV[2];
  71 $maxlong=$ARGV[3];
  72
  73 #$minlat2=int($minlat*2)/2;
  74 #$minlong2=int($minlong*2)/2;
  75 #$maxlat2=ceil($maxlat*2)/2;
  76 #$maxlong2=ceil($maxlong*2)/2;
  77
  78 $minlat2=$minlat;
  79 $minlong2=$minlong;
  80 $maxlat2=$maxlat;
  81 $maxlong2=$maxlong;
  82
  83
  84 #$rect_length=2*abs($maxlong-$minlong) + 2*abs($maxlat-$minlat);
  85 #
  86 #push (@corner_rect2lin, 0);
  87 #push (@corner_rect2lin, abs($maxlong-$minlong));
  88 #push (@corner_rect2lin, abs($maxlong-$minlong) + abs($maxlat-$minlat));
  89 #push (@corner_rect2lin, $rect_length - abs($maxlat-$minlat));
  90
  91 my $osm_id=-1;
  92
  93 my (%new_ways, %nodes, %ways, %tiles, @new_nodes);
  94
  95 my ($i, $j, $way, $node, $tile);
  96
  97 my $debug=0;
  98
  99 sub debug {
 100         print shift if ( $debug );
 101 }
 102
 103 # project a node on a bounding box
 104 sub get_projection {
 105         my ($lat, $lon, $minlat, $minlon, $maxlat, $maxlon) = @_;
 106
 107         my ($a, $b, $c, $d);
 108
 109         $a = abs( $lat - $minlat);
 110         $b = abs( $lat - $maxlat);
 111         $c = abs( $lon - $minlon);
 112         $d = abs( $lon - $maxlon);
 113
 114
 115         if ( $a <= $b ) {
 116                 if ( $a <= $c ) {
 117                         if ( $a <= $d ) {
 118                                 # a
 119                                 $lat=$minlat;
 120                                 #$lon=$node->{"lon"};
 121                         } else {
 122                                 # d
 123                                 #$lat=$node->{"lat"};
 124                                 $lon=$maxlon;
 125                         }
 126                 } else {
 127                         if ( $c <= $d ) {
 128                                 # c
 129                                 #$lat=$node->{"lat"};
 130                                 $lon=$minlon;
 131                         } else {
 132                                 #d
 133                                 #$lat=$node->{"lat"};
 134                                 $lon=$maxlon;
 135                         }
 136                 }
 137         } else {
 138                 if ( $b <= $c ) {
 139                         if ( $b <= $d ) {
 140                                 #b
 141                                 $lat=$maxlat;
 142                                 #$lon=$node->{"lon"};
 143                         } else {
 144                                 # d
 145                                 #$lat=$node->{"lat"};
 146                                 $lon=$maxlon;
 147                         }
 148                 }  else {
 149                         if ( $c <= $d ) {
 150                                 # c
 151                                 #$lat=$node->{"lat"};
 152                                 $lon=$minlon;
 153                         } else {
 154                                 #d
 155                                 #$lat=$node->{"lat};
 156                                 $lon=$maxlon;
 157                         }
 158                 }
 159         }
 160
 161         return (($lat, $lon));
 162 }
 163
 164 # Projet a node on a rectangle to a segment
 165 sub get_proj_rect2lin {
 166         my ($lat, $lon, $minlat, $minlon, $maxlat, $maxlon) = @_;
 167
 168
 169         if ($lat eq $maxlat) {
 170                 return abs($minlon - $lon);
 171         } else {
 172                 if ($lon eq $maxlon) {
 173                         return abs($maxlon-$minlon) + abs($maxlat - $lat);
 174                 } else {
 175                         if ($lat eq $minlat) {
 176                                 return abs($maxlon-$minlon) + abs($maxlat-$minlat) + abs($maxlon-$lon);
 177                         } else {
 178                                 return 2*abs($maxlon-$minlon) + abs($maxlat-$minlat) + abs($minlat - $lat);
 179                         }
 180                 }
 181         }
 182 }
 183
 184 # projet a node from a segment to a rectangle
 185 sub getnode_proj_lin2rect {
 186         my ($length, $minlat, $minlon, $maxlat, $maxlon) = @_;
 187         my ($lat, $lon);
 188
 189         if ( $length > 2*abs($maxlon-$minlon) + abs($maxlat-$minlat)) {
 190                 $lat = ($minlat + $length - 2*abs($maxlon-$minlon) - abs($maxlat-$minlat));
 191                 $lon = ($minlon);
 192         } else {
 193                 if ( $length > abs($maxlon-$minlon) + abs($maxlat-$minlat)) {
 194                         $lat = ($minlat);
 195                         $lon = ($maxlon - $length + abs($maxlon-$minlon) + abs($maxlat-$minlat));
 196                 } else {
 197                         if ( $length > abs($maxlon-$minlon)) {
 198                                 $lat = ($maxlat - $length + abs($maxlon-$minlon));
 199                                 $lon = ($maxlon);
 200                         } else {
 201                                 $lat = ($maxlat);
 202                                 $lon = ($minlon + $length);
 203                         }
 204                 }
 205         }
 206
 207         return (($lat, $lon));
 208 }
 209
 210
 211
 212 # get the distance between to nodes on the same rectangle
 213 sub get_path_length {
 214         my ($lat1, $lon1, $lat2, $lon2, @box) = @_;
 215
 216         my ($a, $b, $rect_length, $dist);
 217
 218
 219         $a = get_proj_rect2lin ( $lat1, $lon1, @box);
 220         $b = get_proj_rect2lin ( $lat2, $lon2, @box);
 221
 222         $dist = $b - $a;
 223         $rect_length = 2*abs( $box[2] - $box[0]) + 2*abs( $box[3] - $box[1]);
 224         return $dist > 0 ? $dist : $dist + $rect_length;
 225 }
 226
 227 # get the nodes to go through to join an other node over a rectangle
 228 sub get_path {
 229         my ($lat1, $lon1, $lat2, $lon2, @box) = @_;
 230
 231         my ($a, $b, $rect_length, $dist, @corner_rect2lin);
 232         my @path=();
 233
 234         $a = get_proj_rect2lin ( $lat1, $lon1, @box);
 235         $b = get_proj_rect2lin ( $lat2, $lon2, @box);
 236
 237         $rect_length = 2*abs( $box[2] - $box[0]) + 2*abs( $box[3] - $box[1]);
 238
 239         $dist = $b - $a;
 240         if ($dist <0) {
 241                 $b+=$rect_length;
 242         }
 243
 244         @corner_rect2lin = ( 0,
 245                 abs( $box[3] - $box[1]),
 246                 abs( $box[3] - $box[1]) + abs( $box[2] - $box[0]),
 247                 $rect_length - abs( $box[2] - $box[0]));
 248
 249         foreach (@corner_rect2lin) {
 250                 if ($_ > $a && $_ < $b) {
 251                         my ($lat, $lon) = getnode_proj_lin2rect ($_, @box);
 252                         push (@path, {"lat" => $lat, "lon" => $lon});
 253                 }
 254         }
 255
 256         foreach (@corner_rect2lin) {
 257                 if ($_ + $rect_length > $a && $_ + $rect_length < $b) {
 258                         my ($lat, $lon) = getnode_proj_lin2rect ($_, @box);
 259                         push (@path, {"lat" => $lat, "lon" => $lon});
 260                 }
 261         }
 262
 263
 264         return @path;
 265 }
 266
 267 # type of tile at zoom12 (land, sea, both, unknown.
 268 # Uses oceantiles_12.dat for tiles@home project
 269 # it checks the tile at level 12 containing the point ($lat, $lon)
 270 # dx and dy are used to specify an offset: read tile (x+dx, y+dy)
 271 # A more clever way whould be to check all the tile contained in and area
 272 # and guess the final type
 273 # one sea or more => sea
 274 # one land or more => land
 275 # only both or unknown => since tiles are quite big, should be land
 276 #
 277 sub lookup_handler {
 278         my ($lat, $lon, $dx, $dy) = @_;
 279
 280         # http://wiki.openstreetmap.org/index.php/Slippy_map_tilenames
 281         # http://almien.co.uk/OSM/Tools/Coord
 282         my $zoom=12;
 283         # $tilex=int( ( $minlon +180)/360 *2**$zoom ) ;
 284         my $tilex=int( ( $lon  +180)/360 *2**$zoom ) ;
 285         # $tiley=int( (1 - log(tan( $maxlat  *pi/180) + sec( $maxlat  *pi/180))/pi)/2 *2**$zoom ) ;
 286         my $tiley=int( (1 - log(tan( $lat  *pi/180) + sec(  $lat  *pi/180))/pi)/2 *2**$zoom ) ;
 287
 288
 289         my $tileoffset = (($tiley+$dy) * (2**12)) + ($tilex+$dx);
 290
 291         my $fh;
 292         open($fh, "<", "data/oceantiles_12.dat") or die ("Cannot open oceantiles_12.dat: $!");
 293
 294
 295         seek $fh, int($tileoffset / 4), 0;
 296         my $buffer;
 297         read $fh, $buffer, 1;
 298         $buffer = substr( $buffer."\0", 0, 1 );
 299         $buffer = unpack "B*", $buffer;
 300         my $str = substr( $buffer, 2*($tileoffset % 4), 2 );
 301
 302         close($fh);
 303
 304         #print "S" if $str eq "10";
 305         #print "L" if $str eq "01";
 306         #print "B" if $str eq "11";
 307         #print "?" if $str eq "00";
 308
 309         return $str;
 310 }
 311
 312
 313 # return the intersection of a segment with the border of a tile
 314 sub intersec_way_tiles {
 315         my ($node_in_lat, $node_in_lon,
 316                 $node_ext_lat,$node_ext_lon,
 317                 $tile_minlat, $tile_minlon, $tile_maxlat, $tile_maxlon) = @_;
 318         my ($lat, $lon);
 319
 320         debug "intersec_way_tiles - $node_in_lat, $node_in_lon && $node_ext_lat,$node_ext_lon && $tile_minlat, $tile_minlon, $tile_maxlat, $tile_maxlon\n";
 321
 322         # case #1: node_ext above tile - try intersec with the top of the tile
 323         if ($node_ext_lat > $tile_maxlat || $node_ext_lat < $tile_minlat) {
 324                 $lat = ($node_ext_lat > $tile_maxlat) ? $tile_maxlat : $tile_minlat;
 325                 if ($node_in_lat != $node_ext_lat) {
 326                         $lon = $node_in_lon + ($lat - $node_in_lat) * ($node_ext_lon - $node_in_lon) / ($node_ext_lat - $node_in_lat);
 327                 } else { # node_ext == node
 328                         $lon = $node_in_lon;
 329                 }
 330                 if ($lon > $tile_minlon && $lon < $tile_maxlon) {
 331                         return ($lat, $lon);
 332                 }
 333         }
 334         # case #2: node ext on the left or on the right of the tile
 335         if ($node_ext_lon > $tile_maxlon || $node_ext_lon < $tile_minlon) {
 336                 $lon = ($node_ext_lon > $tile_maxlon) ? $tile_maxlon : $tile_minlon;
 337                 if ($node_in_lon != $node_ext_lon) {
 338                         $lat = $node_in_lat + ($lon - $node_in_lon) * ($node_ext_lat - $node_in_lat) / ($node_ext_lon - $node_in_lon); # thales
 339                 } else { # node_ext == node
 340                         $lat = $node_in_lat;
 341                 }
 342                 if ($lat > $tile_minlat && $lat < $tile_maxlat) {
 343                         return ($lat, $lon);
 344                 }
 345         }
 346
 347         return ($node_ext_lat, $node_ext_lon);
 348 }
 349
 350 # tells if a node is inside a polygon or not
 351 # not efficient a 100%, but quite good results
 352 # widely inspired from close-areas.pl from the T@h project
 353 sub polygon_contains_node {
 354         my ($way, $node) = @_;
 355
 356
 357         my $counter = 0;
 358
 359         my @way_nodes = $way->nodes();
 360         my $p1 =  $nodes{$way->first_node()};
 361         my $n = scalar(@way_nodes);
 362
 363         for (my $i=1; $i < $n; $i++) {
 364                 my $p2 = $nodes{$way_nodes[$i]};
 365
 366                 if ($node->{"lat"} > $p1->{"lat"} || $node->{"lat"} > $p2->{"lat"}) {
 367                         if ($node->{"lat"} < $p1->{"lat"} || $node->{"lat"} < $p2->{"lat"}) {
 368                                 if ($node->{"lon"} < $p1->{"lon"} || $node->{"lon"} < $p2->{"lon"}) {
 369                                         if ($p1->{"lat"} != $p2->{"lat"}) {
 370                                                 my $xint = ($node->{"lat"}-$p1->{"lat"})*($p2->{"lon"}-$p1->{"lon"})/($p2->{"lat"}-$p1->{"lat"})+$p1->{"lon"};
 371                                                 if ($p1->{"lon"} == $p2->{"lon"} || $node->{"lon"} <= $xint) {
 372                                                         $counter++;
 373                                                 }
 374                                         }
 375                                 }
 376                         }
 377                 }
 378                 $p1 = $p2;
 379         }
 380
 381         return ($counter%2);
 382 }
 383
 384 #################################
 385 print "initialising tiles\n";
 386
 387
 388 debug("there will be ".ceil(abs($maxlong2-$minlong2)/0.5)."x".ceil(abs($maxlat2-$minlat2)/0.5)." tiles\n");
 389
 390 my $nb_tiles_x = ceil(abs($maxlat2-$minlat2)/0.5);
 391 my $nb_tiles_y = ceil(abs($maxlong2-$minlong2)/0.5);
 392
 393 for ($j=0; $j<$nb_tiles_x; $j++) {
 394         for ($i=0; $i<$nb_tiles_y; $i++) {
 395
 396
 397                 my $type = lookup_handler ( $maxlat2-0.5*$j, $minlong2+0.5*$i, 0, 0);
 398
 399                 print "S" if $type eq "10";
 400                 print "L" if $type eq "01";
 401                 print "B" if $type eq "11";
 402                 print "?" if $type eq "00";
 403                 $tiles{$nb_tiles_x * $i+$j} = { "lon" => $minlong2+0.5*$i, "lat" => $maxlat2-0.5*$j, "type" => $type, "id" => ($nb_tiles_x * $i+$j), "crossed" => 0};
 404         }
 405         print "\n";
 406 }
 407
 408 print "\n";
 409
 410
 411
 412 ################################
 413 # find coastlines and store the nodes id
 414 print  "Reading osm source file\n";
 415
 416 my $input;
 417 my $output;
 418 open($input, "<$source") or die "Can't open $source: $!";
 419
 420 my $way_id=0;
 421 my $is_coastline=0;
 422 my @nodes;
 423 while (<$input>) {
 424
 425         if ( /<way id=["'](\d+?)["']/) {
 426                 $way_id=$1;
 427                 $is_coastline=0;
 428                 @nodes=();
 429         }
 430         if ($way_id != 0) {
 431                 if ( /<nd ref=["'](\d+?)["']/) {
 432                         push(@nodes, $1);
 433                 }
 434                 if ( /<tag k=["']natural["'] v=["']coastline["']/) {
 435                         $is_coastline=1;
 436                 }
 437                 if ( /<\/way/ ) {
 438                         # forget node with only one element -> disturb the algo
 439                         if ( $is_coastline && scalar(@nodes) != 1 ) {
 440                                 my $way=way->new($way_id);
 441                                 $way->nodes(@nodes);
 442                                 debug( "found coastline ".$way->id()." from ".$way->first_node()." to ".$way->last_node()." with ".$way->size()." nodes\n");
 443                                 $ways{$way->id()}=$way;
 444                         }
 445
 446                         $way_id=0;
 447                 }
 448         }
 449 }
 450
 451 close($input);
 452
 453
 454
 455
 456 ###########################
 457 # to save memory, only retrieve usefull nodes
 458 print "Retrieving nodes:\n";
 459
 460 $i=0;
 461 foreach (values %ways) {
 462         foreach ($_->nodes()) {
 463                 $nodes{$_}={};
 464                 $i++;
 465         }
 466 }
 467 debug("Created $i nodes\n");
 468
 469 open($input, "<$source") or die "Can't open $source: $!";
 470 open($output, ">$destination") or die "Can't open $destination: $!";
 471
 472 $i=0;
 473 while (<$input>) {
 474         if ( /<node id=['"](\d+?)['"].* lat=['"](.+?)['"] lon=['"](.+?)['"]/ ) {
 475                 if ( exists($nodes{$1}) ) {
 476                         $nodes{$1}->{"id"}=$1;
 477                         $nodes{$1}->{"lat"}=$2;
 478                         $nodes{$1}->{"lon"}=$3;
 479                         $i++;
 480                 }
 481         }
 482         print $output $_ unless ( /<\/osm>/ );
 483 }
 484 close($input);
 485
 486 debug("Retrieved $i nodes\n");
 487
 488
 489
 490 ############################
 491 # If the bounding box given as arguments to the program is smaller than the one used to generate
 492 # the input file, some nodes may be outside -> check for this before adding them to a tile
 493 print "Spliting data into tiles:\n";
 494
 495
 496 foreach $way (values %ways) {
 497         #print "following new way\n";
 498
 499         my $new_way = way->new($osm_id--);
 500         my @way_nodes = ();
 501         my ($tile_id, $new_tile_id);
 502         $tile_id = -1;
 503         my $prev_node;
 504
 505         foreach $node ($way->nodes()) {
 506
 507                 # node out of bbox -> skip the nodes until we reenter the bbox
 508                 if ($nodes{$node}->{"lat"} < $minlat2 || $nodes{$node}->{"lat"} > $maxlat2
 509                         || $nodes{$node}->{"lon"} < $minlong2 || $nodes{$node}->{"lon"} > $maxlong2) {
 510
 511                         # previous node already skiped or 1st node out of the bbox
 512                         if ( $tile_id == -1 ) {
 513                                 $prev_node = $node;
 514
 515
 516                         } else { # previous node in bbox, we add a node at the intersection of the current tile border and the segment [$node, $prev_node]
 517
 518                                 # find the intersection of the current way with the tiles
 519                                 my $tile = $tiles{$tile_id};
 520                                 my ($lat, $lon) = intersec_way_tiles (
 521                                         $nodes{$prev_node}->{"lat"}, $nodes{$prev_node}->{"lon"},
 522                                         $nodes{$node}->{"lat"}, $nodes{$node}->{"lon"},
 523                                         $tile->{"lat"}-0.5, $tile->{"lon"}, $tile->{"lat"}, $tile->{"lon"}+0.5);
 524                                 my $mid_node = {"lat" => $lat, "lon" => $lon, "id" => $osm_id-- };
 525
 526                                 $nodes{ $mid_node->{"id"} } = $mid_node;
 527                                 push( @new_nodes, $mid_node );
 528
 529                                 push ( @way_nodes, $mid_node->{"id"} );
 530
 531                                 # end way here
 532                                 $new_way->nodes(@way_nodes);
 533                                 $new_ways{$new_way->id()} = $new_way;
 534                                 $ways{$new_way->id()} = $new_way;
 535                                 push( @{$tiles{$tile_id}->{"ways"}}, $new_way->id());
 536                                 $tiles{$tile_id}->{"crossed"}=1;
 537
 538                                 # begin new way
 539                                 $new_way = way->new($osm_id--);
 540                                 @way_nodes = ();
 541
 542                                 # reset title_id
 543                                 $tile_id = -1;
 544
 545                         }
 546                 } else {
 547
 548                         my $tile_lat = int( ($maxlat2 - $nodes{$node}->{"lat"} ) / 0.5);
 549                         my $tile_lon = int( ($nodes{$node}->{"lon"} - $minlong2) / 0.5);
 550                         $new_tile_id = $nb_tiles_x * $tile_lon + $tile_lat;
 551
 552                         if ( $tile_id == -1 ) {
 553                                 $tile_id = $new_tile_id;
 554
 555                                 # we have reentered a tile -> add a node on its border at the intersec of the segment [$prev_node, $node]
 556                                 if (defined $prev_node) {
 557                                         # find the intersection of the current way with the tiles
 558                                         my $tile = $tiles{$tile_id};
 559                                         my ($lat, $lon) = intersec_way_tiles (
 560                                                 $nodes{$node}->{"lat"}, $nodes{$node}->{"lon"},
 561                                                 $nodes{$prev_node}->{"lat"}, $nodes{$prev_node}->{"lon"},
 562                                                 $tile->{"lat"}-0.5, $tile->{"lon"}, $tile->{"lat"}, $tile->{"lon"}+0.5);
 563                                         my $mid_node = {"lat" => $lat, "lon" => $lon, "id" => $osm_id-- };
 564
 565                                         $nodes{ $mid_node->{"id"} } = $mid_node;
 566                                         push( @new_nodes, $mid_node );
 567
 568                                         push ( @way_nodes, $mid_node->{"id"} );
 569
 570                                         $tiles{$tile_id}->{"crossed"}=1;
 571                                 }
 572                         }
 573                         # we have reach an other tile
 574                         if ( $new_tile_id != $tile_id ) {
 575
 576                                 my $tile = $tiles{$tile_id};
 577                                 my $new_tile = $tiles{$new_tile_id};
 578
 579                                 # find the intersection of the way with the tiles
 580                                 my ($lat, $lon) = intersec_way_tiles (
 581                                         $nodes{$prev_node}->{"lat"}, $nodes{$prev_node}->{"lon"},
 582                                         $nodes{$node}->{"lat"}, $nodes{$node}->{"lon"},
 583                                         $tile->{"lat"}-0.5, $tile->{"lon"}, $tile->{"lat"}, $tile->{"lon"}+0.5);
 584                                 my $mid_node = {"lat" => $lat, "lon" => $lon, "id" => $osm_id-- };
 585
 586                                 $nodes{ $mid_node->{"id"} } = $mid_node;
 587                                 push( @new_nodes, $mid_node );
 588
 589                                 push ( @way_nodes, $mid_node->{"id"} );
 590
 591                                 $new_way->nodes(@way_nodes);
 592                                 $new_ways{$new_way->id()} = $new_way;
 593                                 $ways{$new_way->id()} = $new_way;
 594                                 push( @{$tiles{$tile_id}->{"ways"}}, $new_way->id());
 595
 596                                 # set both new and old tile as crossed by a way
 597                                 # the new one is important because when we rich the end
 598                                 # of the coast, there will be no "next way" to set it as
 599                                 # crossed
 600                                 $tiles{$tile_id}->{"crossed"}=1;
 601                                 $tiles{$new_tile_id}->{"crossed"}=1;
 602
 603                                 $new_way = way->new($osm_id--);
 604
 605                                 @way_nodes=();
 606
 607
 608                                 # first node on tile border
 609                                 ($lat, $lon) = intersec_way_tiles (
 610                                         $nodes{$node}->{"lat"}, $nodes{$node}->{"lon"},
 611                                         $nodes{$prev_node}->{"lat"}, $nodes{$prev_node}->{"lon"},
 612                                         $new_tile->{"lat"}-0.5, $new_tile->{"lon"}, $new_tile->{"lat"}, $new_tile->{"lon"}+0.5);
 613                                 my $mid_node2 = {"lat" => $lat, "lon" => $lon, "id" => $osm_id-- };
 614
 615                                 $nodes{ $mid_node2->{"id"} } = $mid_node2;
 616                                 push( @new_nodes, $mid_node2 );
 617
 618                                 push ( @way_nodes, $mid_node2->{"id"} );
 619
 620                                 $tile_id = $new_tile_id;
 621                         }
 622
 623                         push (@way_nodes, $node);
 624                         #print "Node $node added to tile $tile_id\n";
 625                         $prev_node=$node;
 626                 }
 627         }
 628
 629         # if new way not empty, add it
 630         if ( scalar(@way_nodes) > 1 && $tile_id != -1) {
 631                 $new_way->nodes(@way_nodes);
 632                 $new_ways{$new_way->id()} = $new_way;
 633                 $ways{$new_way->id()} = $new_way;
 634                 push( @{$tiles{$tile_id}->{"ways"}}, $new_way->id());
 635         }
 636
 637         $way->distroy();
 638
 639 }
 640
 641
 642
 643
 644
 645
 646
 647 ##############################################
 648 print "\n\nWorking on tiles:\n";
 649
 650
 651 #$tile = $tiles{54};
 652 foreach $tile (values %tiles) {
 653
 654         print( "->tile: ".$tile->{"id"}."\n");
 655
 656         my @box = ($tile->{"lat"}-0.5, $tile->{"lon"}, $tile->{"lat"}, $tile->{"lon"}+0.5);
 657
 658         # No way going through the tile -> guess if it is sea or lanf
 659         # SHould check if it containds nodes -> type unknown + nodes should mean sea + islands
 660         unless ($tile->{"crossed"}) {
 661                 my $tile_is_sea=0;
 662
 663                 #  00 => unknown
 664                 #  01 => land
 665                 #  10 => sea
 666                 #  11 => coast
 667
 668                 # water
 669                 if ($tile->{"type"} eq "10") {
 670                         $tile_is_sea=1;
 671
 672                 } elsif ($tile->{"type"} ne "01") { # not land
 673
 674                         my %temp = ("00"=>0, "10"=>0, "01"=>0, "11"=>0);
 675                         $temp{lookup_handler($tile->{"lat"}, $tile->{"lon"}, -1, 0)}++;
 676                         $temp{lookup_handler($tile->{"lat"}, $tile->{"lon"}, +1, 0)}++;
 677                         $temp{lookup_handler($tile->{"lat"}, $tile->{"lon"}, 0, -1)}++;
 678                         $temp{lookup_handler($tile->{"lat"}, $tile->{"lon"}, 0, +1)}++;
 679
 680                         if( $temp{"10"} > $temp{"01"} ) {
 681                                 $tile_is_sea=1;
 682                         } elsif ( ($tile->{"type"} eq "11") and ( $temp{"01"} == 0 ) ) {
 683                                 # if the tile is marked coast but no land near, assume it's a group of islands instead of lakes.
 684                                 $tile_is_sea=1;
 685                         } # else = land
 686
 687                 }
 688
 689                 # add nodes and a way to draw the sea all over the tile
 690                 if ($tile_is_sea) {
 691                         my @way_nodes_id = ( $osm_id, $osm_id-1, $osm_id-2, $osm_id-3);
 692                         my @way_nodes = (
 693                                 { "lat" => $tile->{"lat"}, "lon" => $tile->{"lon"}, "id" => $osm_id--},
 694                                 { "lat" => $tile->{"lat"}, "lon" => $tile->{"lon"}+0.5, "id" => $osm_id--},
 695                                 { "lat" => $tile->{"lat"}-0.5, "lon" => $tile->{"lon"}+0.5, "id" => $osm_id--},
 696                                 { "lat" => $tile->{"lat"}-0.5, "lon" => $tile->{"lon"}, "id" => $osm_id--});
 697
 698                         my $way = way->new($osm_id--);
 699                         $way->nodes(@way_nodes_id);
 700
 701                         push( @new_nodes, @way_nodes );
 702                         foreach (@way_nodes) {
 703                                 $nodes{$_->{"id"}} = $_;
 704                         }
 705
 706                         $new_ways{$way->id()} = $way;
 707                         $ways{$way->id()} = $way;
 708                         push( @{$tile->{"ways"}}, $way->id());
 709                 }
 710         }
 711
 712         # next tile if there is nothing in this one
 713         next if ( ! $tile->{"ways"} );
 714
 715         #========================================
 716         # merge ways that can be merged
 717         debug "\tBuilding tree:\n";
 718
 719         foreach (@{$tile->{"ways"}}) {
 720                 $way = $new_ways{$_};
 721                 debug( "node ".$way->id()." (size: ".$way->size().")\n" );
 722                 # only for not yet complete polygons and ways without prev and next
 723                 if ($way->next() == 0 || $way->prev() == 0) {
 724                         foreach (@{$tile->{"ways"}}) {
 725                                 my $way2=$new_ways{$_};
 726                                 if ($way->last_node() == $way2->first_node() ) {
 727                                         if ($way->next() == 0) {
 728                                                 debug( "\tnext = ".$way2->id()."\n");
 729                                                 $way->next($way2->id());
 730                                         } else {
 731                                                 debug( "\t error: next way duplicated: old:".$way->next()." new:".$way2->id()." - node: ".$way->last_node()."\n");
 732                                         }
 733                                 }
 734                                 if ($way->first_node() == $way2->last_node() ) {
 735                                         if ($way->prev() == 0) {
 736                                                 debug( "\tprev = ".$way2->id()."\n");
 737                                                 $way->prev($way2->id());
 738                                         } else {
 739                                                 debug( "\t error: prev way duplicated: old:".$way->prev()." new:".$way2->id()." - node: ".$way->first_node()."\n");
 740                                         }
 741                                 }
 742                         }
 743                 } else {
 744                         debug( "\tAlready ok\n" );
 745                 }
 746         }
 747         #=============================
 748         # if some ways are not closed, we need to close them.
 749         # It is done by adding new nodes on the border of the tile
 750         # and joining them
 751         debug ("\n\nClosing loops:\n");
 752
 753         my (@first_nodes, @last_nodes);
 754
 755         # project nodes on the border of the tile
 756         # Some nodes are already on the border of the tile (mainly because of
 757         # intersec_way_tiles)
 758         foreach (@{$tile->{"ways"}}) {
 759                 my $way = $new_ways{$_};
 760
 761                 next if ( $way->first_node() == $way->last_node()); # already a loop
 762
 763                 if ( $way->prev() == 0) {
 764                         my @way_nodes = $way->nodes();
 765
 766                         my $lat = $nodes{$way_nodes[0]}->{"lat"};
 767                         my $lon = $nodes{$way_nodes[0]}->{"lon"};
 768                         ($lat, $lon) = get_projection( $lat, $lon, $tile->{"lat"}-0.5, $tile->{"lon"}, $tile->{"lat"}, $tile->{"lon"}+0.5 );
 769                         my $node = { "lat" => $lat, "lon" => $lon, "id" => ($osm_id--)};
 770                         unshift( @way_nodes, $node->{"id"} );
 771                         push( @new_nodes, $node );
 772                         $nodes{$node->{"id"}}=$node;
 773                         push( @first_nodes, $node->{"id"});
 774
 775                         $way->nodes(@way_nodes);
 776                 }
 777                 if ( $way->next() == 0) {
 778
 779                         my @way_nodes = $way->nodes();
 780
 781                         my $lat = $nodes{$way_nodes[ @way_nodes - 1 ]}->{"lat"};
 782                         my $lon = $nodes{$way_nodes[ @way_nodes - 1 ]}->{"lon"};
 783                         ($lat, $lon) = get_projection( $lat, $lon, $tile->{"lat"}-0.5, $tile->{"lon"}, $tile->{"lat"}, $tile->{"lon"}+0.5 );
 784                         my $node = { "lat" => $lat, "lon" => $lon, "id" => $osm_id--};
 785                         push( @way_nodes, $node->{"id"} );
 786                         push( @new_nodes, $node );
 787                         $nodes{$node->{"id"}}=$node;
 788                         push( @last_nodes, $node->{"id"});
 789
 790                         $way->nodes(@way_nodes);
 791
 792                 }
 793
 794         }
 795
 796         # link nodes -> find a way around the tile
 797
 798
 799         my $last;
 800         foreach $last (@last_nodes) {
 801                 debug "Looking for path ...\n";
 802                 my ($path1, $path2);
 803                 $path1=0;
 804
 805                 # find the sortest path between nodes
 806                 # -> will give the next node around the tile
 807                 my $first;
 808                 my $first_elected;
 809                 foreach $first (@first_nodes) {
 810                         $path2=get_path_length($nodes{$last}->{"lat"},
 811                                 $nodes{$last}->{"lon"},
 812                                 $nodes{$first}->{"lat"},
 813                                 $nodes{$first}->{"lon"},
 814                                 @box);
 815                         if ($path1 == 0 || $path1 > $path2) {
 816                                 $path1=$path2;
 817                                 $first_elected=$first;
 818                         }
 819                 }
 820
 821                 # get the path around the tile between the two nodes
 822                 my @path = get_path ( $nodes{$last}->{"lat"},
 823                         $nodes{$last}->{"lon"},
 824                         $nodes{$first_elected}->{"lat"},
 825                         $nodes{$first_elected}->{"lon"},
 826                         @box);
 827
 828                 foreach (@path) {
 829                         $_->{"id"} = $osm_id--;
 830                         push(@new_nodes, $_);
 831                         $nodes{$_->{"id"}} = $_;
 832                 }
 833
 834
 835                 unshift(@path, $nodes{$last});
 836                 push(@path, $nodes{$first_elected});
 837
 838                 my $way = way->new($osm_id--);
 839                 my @way_nodes=();
 840                 debug( "new way ".$way->id().":\n");
 841                 foreach (@path) {
 842                         push(@way_nodes, $_->{"id"});
 843                         debug( "nodes id: ".$_->{"id"}."\tlat: ".$_->{"lat"}." - long: ".$_->{"lon"}."\n");
 844                 }
 845                 $way->nodes(@way_nodes);
 846
 847                 push( @{$tiles{$tile->{"id"}}->{"ways"}}, $way->id());
 848                 $new_ways{$way->id()} = $way;
 849                 $ways{$way->id()} = $way;
 850
 851
 852
 853         }
 854
 855
 856         # linking the ways
 857         # This time all ways should be loops because of the ways we added before
 858         debug "\n\nBuilding tree with new ways:\n";
 859
 860         foreach (@{$tile->{"ways"}}) {
 861                 my $way = $new_ways{$_};
 862                 debug( "node ".$way->id()." (size: ".$way->size().")\n" );
 863                 # only for not yet complete polygons and ways without prev and next
 864                 if ($way->next() == 0 || $way->prev() == 0) {
 865                         foreach (@{$tile->{"ways"}}) {
 866                                 my $way2=$new_ways{$_};
 867                                 if ($way->last_node() == $way2->first_node() ) {
 868                                         if ($way->next() == 0) {
 869                                                 debug( "\tnext = ".$way2->id()."\n");
 870                                                 $way->next($way2->id());
 871                                         } else {
 872                                                 debug( "\t error: next way duplicated: old:".$way->next()." new:".$way2->id()." - node: ".$way->last_node()."\n");
 873                                         }
 874                                 }
 875                                 if ($way->first_node() == $way2->last_node() ) {
 876                                         if ($way->prev() == 0) {
 877                                                 debug( "\tprev = ".$way2->id()."\n");
 878                                                 $way->prev($way2->id());
 879                                         } else {
 880                                                 debug( "\t error: prev way duplicated: old:".$way->prev()." new:".$way2->id()." - node: ".$way->first_node()."\n");
 881                                         }
 882                                 }
 883                         }
 884                 } else {
 885                         debug( "\tAlready ok\n" );
 886                 }
 887         }
 888
 889
 890         #merge ways to have only simple ways that are circular
 891         debug "merge ways\n";
 892
 893         # copy into @a because we add new ways to the array inside the loop
 894         # -> otherwise it creates an infinite loop
 895         my @a = @{$tile->{"ways"}};
 896         foreach (@a) {
 897 #               print "way $_\n";
 898                 my $way = $new_ways{$_};
 899                 next unless ($way->is_alive());
 900
 901                 my @way_nodes;
 902                 while ($way->is_alive()) {
 903
 904                         push (@way_nodes, $way->nodes());
 905
 906                         # Distroy the way -> it won't be counted twince thus this will stop the loop
 907                         # and it will free some memory
 908                         $way->distroy();
 909                         last if ($way->next() == 0);
 910                         $way=$new_ways{ $way->next() };
 911                 }
 912                 my $new_way = way->new($osm_id--);
 913                 $new_way->nodes(@way_nodes);
 914
 915                 push( @{$tiles{$tile->{"id"}}->{"ways"}}, $new_way->id());
 916                 $new_ways{$new_way->id()} = $new_way;
 917                 $ways{$new_way->id()} = $new_way;
 918
 919         }
 920
 921
 922         ## sort the ways between sea and land
 923         foreach (@{$tile->{"ways"}}) {
 924                 my $way = $new_ways{$_};
 925                 next unless ($way->is_alive());
 926 #                               print "---- ".$way->id().": ".$way->first_node()." -> ".$way->last_node()."\n";
 927                 foreach (@{$tile->{"ways"}}) {
 928                         my $way2 = $new_ways{$_};
 929                         # - don't check if a way is included in itself
 930                         # - check only for ways of type "outer" which is the default type
 931                         #    => actualy should ckeck it. if a not is inside a "inner" way then it belongs
 932                         #    to a lake that is in an island
 933                         if ($way2->is_alive() && $way2->id() != $way->id() && $way2->type() eq "outer") {
 934                                 if (polygon_contains_node($way2, $nodes{$way->first_node()})) {
 935                                         $way->type("inner");
 936                                         last;
 937                                 }
 938                         }
 939                 }
 940         }
 941 }
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954 ########################################
 955 print "\n\nGenerating new data:\n";
 956
 957 # add new nodes to output file
 958 foreach $node (@new_nodes) {
 959         debug( "adding node ".$node->{"id"}." to osm file\n");
 960         print $output "  <node id=\"".$node->{"id"}."\" timestamp=\"2008-07-02T22:57:53Z\" user=\"latouche\" lat=\"".$node->{"lat"}."\" lon=\"".$node->{"lon"}."\"/>\n";
 961 }
 962
 963 # add new ways (highway=motorway for debug)
 964 #foreach $way (values %new_ways) {
 965 #       print $output "  <way id=\"".$way->id()."\" timestamp=\"2008-07-02T22:57:53Z\" user=\"latouche\">\n";
 966 #       foreach $node ($way->nodes()) {
 967 #               print $output "    <nd ref=\"".$node."\"/>\n" if defined $node;
 968 #       }
 969 #
 970 #       print $output "    <tag k=\"highway\" v=\"motorway\"/>\n";
 971 #       print $output "  </way>\n";
 972 #
 973 #}
 974
 975
 976 # add new ways to output file
 977 foreach $tile (values %tiles) {
 978         print( "Writing tile: ".$tile->{"id"}."\n");
 979         foreach (@{$tile->{"ways"}}) {
 980                 my $way = $new_ways{$_};
 981                 my $type = "inner";
 982                 next unless ($way->is_alive());
 983                 print $output "  <way id=\"".($osm_id--)."\" timestamp=\"2008-07-02T22:57:53Z\" user=\"latouche\">\n";
 984                 while ($way->is_alive()) {
 985
 986                         $type = ($way->type() eq "outer") ? "water" : "land";
 987
 988                         foreach $node ($way->nodes()) {
 989                                 print "undef node in ".$way->id()."\n" unless defined $node;
 990                                 print $output "    <nd ref=\"".$node."\"/>\n" if defined $node;
 991                         }
 992                         $way->distroy();
 993                         last if ($way->next() == 0);
 994                         $way=$new_ways{ $way->next() };
 995                 }
 996                 print $output "    <tag k=\"natural\" v=\"".$type."\"/>\n";
 997                 print $output "  </way>\n";
 998         }
 999 }
1000
1001
1002 print $output "</osm>\n";
1003
1004 close($output);
1005
1006 exit(0);