source: subversion/applications/rendering/tilesAtHome-dev/trunk/close-areas.pl @ 23767

Last change on this file since 23767 was 14988, checked in by frederik, 11 years ago

fix id/uid ambiguity in close-areas node parsing regex

  • Property svn:eol-style set to native
File size: 31.2 KB
Line 
1#!/usr/bin/perl
2
3# A program to form closed areas from incomplete sets of bordering
4# segments, by making an assumption about which side of the segment
5# the enclosed area lies on.
6#
7# Currently only supports "coastline" ways but may be extended in
8# the future.
9#
10# A detailed discussion should be in the Wiki:
11# http://wiki.openstreetmap.org/index.php/Tiles%40home/Dev/Interim_Coastline_Support
12#
13# Written by Frederik Ramm <frederik@remote.org> - public domain
14
15use lib './lib';
16use Math::Trig;
17use Math::Complex;
18use tahproject;
19use strict;
20
21my $nodecount = 0;
22my $waycount = 0;
23my $relcount = 0;
24my $halfpi = pi()/2;
25my $twopi = pi()*2;
26
27my $LimitY = pi();
28my $LimitY2 = -pi();
29my $RangeY = $twopi;
30my $debug=0;
31my $make_open_ways=0;
32my $segments;
33my $nodes;
34
35# The direction in which areas are closed - cw (clockwise), or ccw
36# (counter-clockwise). For OSM we normally use cw, which means that
37# area boundaries are drawn such that the area is "to the right" of
38# the segment (i.e. coastline segments have water on their right).
39#
40# "ccw" is unused but supported in case somebody needs it.
41my $direction="cw";
42
43my @coastline_segments;
44
45my $tilex;
46my $tiley; 
47my $zoom;
48my $minlat;
49my $minlon;
50my $maxlat;
51my $maxlon;
52my $border_crossed;
53
54if (scalar(@ARGV) == 2 or scalar(@ARGV) == 3)
55{
56    ($tilex, $tiley, $zoom) = @ARGV;
57    if (!$zoom) {$zoom = 12;}
58    ($maxlat, $minlat) = Project($tiley, $zoom);
59    ($minlon, $maxlon) = ProjectL($tilex, $zoom);
60}
61elsif (scalar(@ARGV) == 4)
62{
63    ($minlat, $minlon, $maxlat, $maxlon) = @ARGV;
64}
65else
66{
67    die "need either 'tilex tiley [zoom]' (defaults to zoom=12) or 'minlat minlon maxlat maxlon' on cmdline";
68}
69
70my $sigma = 0.01;
71
72my $helpernodes =
73    { 0 => [ $maxlat + $sigma, $maxlon + $sigma], 
74      1 => [ $maxlat + $sigma, $minlon - $sigma ], 
75      2 => [ $minlat - $sigma, $minlon - $sigma ], 
76      3 => [ $minlat - $sigma, $maxlon + $sigma ] };
77
78my $copy = 1;
79my $waybuf;
80my $is_coastline;
81my @seglist;
82my $last_node_ref; 
83my $segcount = 0; 
84
85while(<STDIN>)
86{
87    while(/(<[^'"<>]*((["'])[^\3]*?\3[^<>"']*)*>)/og)
88    {
89        my $xmltag = $1;
90        if ($xmltag =~ /^\s*<node.*\sid=["'](\d+)['"].*lat=["']([0-9.Ee-]+)["'].*lon=["']([0-9.Ee-]+)["']/)
91        {
92            $nodes->{$1} = { "lat" => $2, "lon" => $3 };
93        }
94        elsif ($xmltag =~ /^\s*<way /)
95        {
96            $copy = 0;
97            undef $waybuf;
98            undef @seglist;
99            undef $is_coastline;
100            undef $last_node_ref;
101        }
102        elsif($xmltag =~ /^\s*<\/osm/)
103        {
104            last;
105        }
106        elsif($xmltag =~ /^\s*<(osm.*)>/)
107        {
108            # If frollo encounters an empty file, it outputs <osm foo />. Detect this and exit
109            if (substr($1, -1) eq "/" )
110            { print $xmltag; exit }
111           
112            if (/version=['"](.*?)['"]/)
113            {
114                die ("close-areas.pl does not support version $1") unless ( ($1 eq "0.5") || ($1 eq "0.6") );
115            }
116        }
117
118        if ($copy)
119        {
120            print $xmltag."\n";
121        }
122        elsif ($xmltag =~ /^\s*<tag k=['"]natural["'].*v=["']coastline['"]/)
123        {
124            $is_coastline = 1; 
125        }
126        elsif ($xmltag =~ /^\s*<nd ref=['"](\d+)["']/)
127        {
128            $waybuf .= $xmltag . "\n";
129            if (defined($last_node_ref))
130            {
131                push(@seglist, { "from" => $last_node_ref, "to" => $1 });
132            }
133            $last_node_ref = $1;
134        }
135        elsif ($xmltag =~ /^\s*<\/way/)
136        {
137            $copy = 1;
138            # non-coastal ways are written and forgotten
139            if (!$is_coastline)
140            {
141                print $waybuf;
142                print $xmltag . "\n";
143            }
144            # coastal ways are forgotten too, but their segments are kept.
145            else
146            {
147                foreach my $seg(@seglist)
148                {
149                    $segments->{++$segcount} = $seg;
150                }
151            }
152        }
153        else
154        {
155            $waybuf .= $xmltag . "\n";
156        }
157    }
158}
159
160# file fully read. delete non-coastal segments (we have printed them
161# all, so no need for us to keep them), and check coastal segments for
162# intersection with the bounding box.
163
164foreach my $seg(keys(%$segments))
165{
166    my $fromnode = $nodes->{$segments->{$seg}->{"from"}};
167    my $tonode = $nodes->{$segments->{$seg}->{"to"}};
168
169    if (!(defined($fromnode) && defined($tonode)))
170    {
171        delete $segments->{$seg};
172        print "delete segment - incomplete $seg\n" if ($debug);
173        next;
174    }
175
176    # returns 0, 1, or 2 points.
177    # (may probably return 3 or 4 points in freak cases)
178    my $intersect = compute_bbox_intersections($fromnode, $tonode);
179    printf "intersections for $seg: %d\n", scalar(@$intersect) if ($debug);
180
181    if (!scalar(@$intersect))
182    {
183        # this segment has no intersections with the bounding box.
184        if (!node_is_inside($fromnode))
185        {
186            # this segment is fully outside the bounding box, and of
187            # no concern to us.
188            delete $segments->{$seg};
189            print "delete $seg fully out\n" if ($debug);
190        }
191    }
192    elsif (scalar(@$intersect) == 1)
193    {
194        # this segments enters OR exits the bounding box. find out which,
195        # and tag accordingly.
196        if (node_is_inside($fromnode))
197        {
198            $segments->{$seg}->{"exit_intersect"} = $intersect->[0];
199            $segments->{$seg}->{"exit_intersect_angle"} =
200                compute_angle_from_bbox_center($intersect->[0]);
201        }
202        else
203        {
204            $segments->{$seg}->{"entry_intersect"} = $intersect->[0];
205            $segments->{$seg}->{"entry_intersect_angle"} =
206                compute_angle_from_bbox_center($intersect->[0]);
207        }
208    }
209    else
210    {
211        # this segment enters AND exits the bounding box. as intersection
212        # points are ordered by distance from the segment's origin, we
213        # assume that the very first intersection is the entry point and
214        # the very last is the exit point.
215        #   
216        # FIXME: segments like this - probably a very long one cutting right
217        # through the box or a short one cutting diagonally at one edge -
218        # are very little tested.
219        $segments->{$seg}->{"entry_intersect_angle"} =
220            compute_angle_from_bbox_center($intersect->[0]);
221        $segments->{$seg}->{"exit_intersect_angle"} =
222            compute_angle_from_bbox_center($intersect->[scalar(@$intersect)-1]);
223        $segments->{$seg}->{"entry_intersect"} = $intersect->[0];
224        $segments->{$seg}->{"exit_intersect"} = $intersect->[scalar(@$intersect)-1];
225    }
226}
227
228# if no coastline segments are present, switch over to
229# special handler to decide whether to draw a blue tile.
230
231if (!scalar(%$segments))
232{
233    ## instead of filling the tile here, we do it at the end...
234    $border_crossed = 0;
235    goto ENDOSM;
236}
237
238# now start building artificial ways for coastline segments.
239#
240# strategy:
241# 1. grab any one segment. if none available, we're done.
242# 2. find segment that starts where previous segment ends,
243#    repeat until none found OR back where we started.
244# 3. if back where we started, write the circular list of segments,
245#    delete them, and repeat from 1.
246# 4. if none found AND last position is inside drawing area,
247#    remove processed segments without writing, and repeat from 1.
248# 5. if none found AND last position is outside drawing area,
249#    search clockwise outside of drawing area for another
250#    segment that begins here, create articifial segment joining
251#    latest position and segment's start node, and continue with
252#    step 2.
253#
254# Special rule for creating artificial segments: artificial segments
255# must never intersect with drawing area. use artificial nodes outside
256# the four corners of the drawing area to route segments if necessary.
257
258# copy list of segment ids
259my $available_segs;
260grep { $available_segs->{$_} = 1 } keys(%$segments);
261
262while (scalar(%$available_segs))
263{
264    my @seglist;
265    # grab any segment as the first one
266    my @tmp = keys(%$available_segs);
267    my $currentseg = shift(@tmp);
268    printf "GRAB %d (from %d to %d)\n", $currentseg, $segments->{$currentseg}->{"from"}, $segments->{$currentseg}->{"to"} if ($debug);
269    delete $available_segs->{$currentseg};
270    #printf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
271    my $currentnode = $segments->{$currentseg}->{"to"};
272    push (@seglist, $currentseg);
273
274STRING:
275    # find a segment that begins where the previous one ended
276    while ($currentseg)
277    {
278        printf "SEARCH begin at %d\n", $currentnode if ($debug);
279        undef $currentseg;
280        foreach my $seg(keys(%$available_segs))
281        {
282            #printf "  TEST seg %d begin at %d\n", $seg, $segments->{$seg}->{"from"} if ($debug);
283            if ($segments->{$seg}->{"from"} == $currentnode)
284            {
285                printf "  FOUND seg %d begin at %d\n", $seg, $segments->{$seg}->{"from"} if ($debug);
286                push (@seglist, $seg);
287                $currentseg = $seg;
288                $currentnode = $segments->{$currentseg}->{"to"};
289                delete $available_segs->{$seg};
290                #printf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
291                last;
292            }
293        }
294    }
295
296    # no more segments found. do we have a loop?
297    if ($currentnode == $segments->{$seglist[0]}->{"from"})
298    {
299        printf("LOOP\n") if ($debug);
300        # loop detected. store the segment list for later output,
301        # and move on trying to connect the rest.
302        push(@coastline_segments, \@seglist);
303        next;
304    }
305
306    my $lastseg = @seglist[scalar(@seglist)-1];
307    my $exit_angle = $segments->{$lastseg}->{"exit_intersect_angle"};
308
309    # are we inside the drawable area? (we are if the last segment did not
310    # exit the drawable area.)
311    # if yes, give up as we obviously have an imcomplete piece of coastline.
312    if (!defined($exit_angle))
313    {
314        printf("NOT OUTSIDE\n") if ($debug);
315        # this is a debug option that allows one to detect the incomplete
316        # way by looking at the output file with JOSM etc.; it is not intended
317        # for production use!
318        if ($make_open_ways)
319        {
320            make_way(\@seglist, 1);
321        }
322        next;
323    }
324
325    # "else" case: we are outside the drawable area and want
326    # to find another segment that begins outside the drawable
327    # area and enters it.
328   
329    my $segs_to_check;
330    $segs_to_check = [];
331    foreach my $seg(keys(%$available_segs))
332    {
333        push(@$segs_to_check, $seg) if defined($segments->{$seg}->{"entry_intersect"});
334    }
335
336    # we will also accept the first segment of the current way
337    # if it is outside. this is a special case as, being used
338    # already, it isn't in the $segments list any
339    # more
340   
341    push(@$segs_to_check, $seglist[0]) 
342        if defined($segments->{$seglist[0]}->{"entry_intersect"});
343
344    # sort all segments entering the drawable area by angle from area
345    # centrepoint to the point where they enter the area.
346    my @sorted_segs_to_check;
347    @sorted_segs_to_check = sort
348        { $segments->{$a}->{"entry_intersect_angle"} <=> 
349          $segments->{$b}->{"entry_intersect_angle"} } @$segs_to_check;
350
351    @sorted_segs_to_check = reverse @sorted_segs_to_check if ($direction eq "cw");
352
353    # find the nearest entering segment.
354    my $found;
355    $found = 0;
356    foreach my $seg(@sorted_segs_to_check)
357    {
358        if ($direction eq "cw")
359        {
360            next if ($segments->{$seg}->{"entry_intersect_angle"}) > $exit_angle;
361        }
362        else
363        {
364            next if ($segments->{$seg}->{"entry_intersect_angle"}) < $exit_angle;
365        }
366        printf("use seg %d angle %f\n", $seg, $segments->{$seg}->{"entry_intersect_angle"}) if ($debug);
367        $found = $seg;
368        last;
369    }
370    if (!$found)
371    {
372        foreach my $seg(@sorted_segs_to_check)
373        {
374            printf("use (2) seg %d angle %f\n", $seg, $segments->{$seg}->{"entry_intersect_angle"}) if ($debug);
375            $found = $seg;
376            last;
377        }
378    }
379
380    # if no segment remains outside, give up
381    if (!$found)
382    {
383        printf("NO SEG OUTSIDE\n") if ($debug);
384        # this is a debug option that allows one to detect the incomplete
385        # way by looking at the output file with JOSM etc.; it is not intended
386        # for production use!
387        if ($make_open_ways)
388        {
389            make_way(\@seglist, 1);
390        }
391        next;
392    }
393   
394    # at this point, we have a list of segments that ends with a segment
395    # leaving the drawable area, and we also have the next segment where
396    # the coastline comes back into the drawable area. we need to connect
397    # them with an artifical segment - or more than one.
398
399    # If a coastline leaves the visible area at the top (side 1) and comes
400    # back in from the right (side 0), then we must not make a direct connection
401    # between the points for fear of clipping our viewport; instead, we must
402    # "hop" over the top right corner. Same for other sides. The "helpernodes"
403    # hash contains nodes to be used for hopping from one side to the next.
404    #
405    # (an extreme case of this is coastline leaving and entering at the same
406    # side but leaving south of where it enters - need to go around all four
407    # corners then, for a total of 5 artificial segments.)
408
409    $border_crossed = 1;
410    my $height = $maxlat-$minlat;
411    my $width = $maxlon-$minlon;
412
413    # exit_angle is already defined
414    my $entry_angle = $segments->{$found}->{"entry_intersect_angle"};
415
416    my $exit_side = $segments->{$lastseg}->{"exit_intersect"}->{"side"};
417    my $entry_side = $segments->{$found}->{"entry_intersect"}->{"side"};
418
419    printf("exit angle %s entry angle %s\n", $exit_angle, $entry_angle) if $debug;
420    printf("exit side %d entry side %d\n", $exit_side, $entry_side) if $debug;
421
422    # the following two blocks (similar but not identical for clockwise/
423    # counter-clockwise) will find out whether we need to "go around corners"
424    # and add 0-4 segments if needed.
425    if ($direction eq "ccw")
426    {
427        # $min_once is for the special case where entry and exit side are
428        # identical but we still need to go all the way around the box.
429        my $min_once;
430        if ($exit_side == 0 and $entry_side == 0) 
431        {
432            # Take into account that the angle flips from 2*pi() to zero at this side
433            $min_once = (($exit_angle > $entry_angle and ($exit_angle - $entry_angle) < pi()) ||
434                         ($exit_angle < pi() and $entry_angle > pi()));
435        }
436        else 
437        {
438            $min_once = ($exit_angle > $entry_angle); 
439        }
440        printf("min_once=%d\n", $min_once) if $debug;
441        for (my $i = $exit_side;; $i++)
442        {
443            $i=0 if ($i==4);
444            last if ($i==$entry_side && !$min_once);
445            $min_once = 0;
446            last if (check_segments_intersect($helpernodes->{$i}->[1], $helpernodes->{$i}->[0], $nodes->{$currentnode}->{"lon"}, $nodes->{$currentnode}->{"lat"}, $segments) );
447            my $newnode = make_node($helpernodes->{$i}->[0], $helpernodes->{$i}->[1]);
448            my $newseg = make_seg($currentnode, $newnode);
449            $currentnode = $newnode;
450            push(@seglist, $newseg);
451        }
452    }
453    else
454    {
455        my $min_once;
456        if ($exit_side == 0 and $entry_side == 0)
457        {
458            $min_once = (($exit_angle < $entry_angle and ($entry_angle - $exit_angle) < pi()) ||
459                         ($exit_angle > pi() and $entry_angle < pi()));
460        }
461        else
462        {
463            $min_once = ($exit_angle < $entry_angle);
464        }
465        printf("min_once=%d\n", $min_once) if $debug;
466        for (my $i = $exit_side;; $i--)
467        {
468            $i=3 if ($i==-1);
469            last if ($i==$entry_side && !$min_once);
470            $min_once = 0;
471            my $helper = $i-1;
472            $helper = 3 if ($helper == -1);
473            last if (check_segments_intersect($helpernodes->{$helper}->[1], $helpernodes->{$helper}->[0], $nodes->{$currentnode}->{"lon"}, $nodes->{$currentnode}->{"lat"}, $segments) );
474            my $newnode = make_node($helpernodes->{$helper}->[0], $helpernodes->{$helper}->[1]);
475            my $newseg = make_seg($currentnode, $newnode);
476            $currentnode = $newnode;
477            push(@seglist, $newseg);
478        }
479    }
480
481    my $newseg = make_seg($currentnode, $segments->{$found}->{"from"});
482    push(@seglist, $newseg);
483
484    # if the segment we have found is our start segment (which we added
485    # as a special case earlier!), we have a closed way.
486    if ($found == $seglist[0])
487    {
488        printf("CLOSED\n") if ($debug);
489        push(@coastline_segments, \@seglist);
490        next;
491    }
492
493    # else we just plod on.
494    push (@seglist, $found);
495    $currentseg = $found;
496    $currentnode = $segments->{$found}->{"to"};
497    delete $available_segs->{$found};
498    #nprintf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
499    goto STRING;
500}
501
502ENDOSM:
503
504# if we had no coastline intersection with the bounding box, but coastline
505# was present, we have an island situation and need to add a blue background.
506# FIXME what if we have a "little lake" situation?
507
508unless( $border_crossed )
509{
510  my $state = lookup_handler($helpernodes, $tilex, $tiley, $zoom);
511  if( $state eq "10" )
512  {
513    # sea
514    addBlueRectangle($helpernodes);
515  }
516  elsif ( $state eq "01" )
517  {
518    #land
519  }
520  else # state 00 or 11 (unknown or mixed)
521  { 
522    my %temp = ("00"=>0, "10"=>0, "01"=>0, "11"=>0);;
523    $temp{lookup_handler($helpernodes, $tilex-1, $tiley, $zoom)}++;
524    $temp{lookup_handler($helpernodes, $tilex+1, $tiley, $zoom)}++;
525    $temp{lookup_handler($helpernodes, $tilex, $tiley-1, $zoom)}++;
526    $temp{lookup_handler($helpernodes, $tilex, $tiley+1, $zoom)}++;
527
528    if( $temp{"10"} + $temp{"11"} > $temp{"01"} ) # if more sea/coast than land around.
529    {
530      addBlueRectangle($helpernodes);
531    }
532    elsif ( ($state eq "11") and ( $temp{"01"} == 0 ) ) 
533    # if the tile is marked coast but no land near, assume it's a group of islands instead of lakes.
534    {
535      # coast
536      addBlueRectangle($helpernodes);
537    }
538    else
539    {
540      #land
541    }
542  }
543}
544
545if (scalar(@coastline_segments) == 1)
546{
547    make_way($coastline_segments[0]);
548}
549else
550{
551    make_multipolygon(\@coastline_segments);
552}
553print "</osm>\n";
554
555
556
557
558# checks segment from ($x0,$y0) to ($x1,$y1) for intersection with all segments
559# in $s
560sub check_segments_intersect
561{
562        my ($x0, $y0, $x1, $y1, $s)= @_;
563       
564        printf("check_segments_intersect: (%f,%f)  (%f,%f)   %s\n", $x0, $y0, $x1, $y1, join(", ",keys(%$s)) ) if ($debug);
565        foreach my $seg(keys(%$s))
566        {
567                my $from= $s->{$seg}->{"from"};
568                my $to  = $s->{$seg}->{"to"};
569                my $x2= $nodes->{$from}->{"lon"};
570                my $y2= $nodes->{$from}->{"lat"};
571                my $x3= $nodes->{$to}  ->{"lon"};
572                my $y3= $nodes->{$to}  ->{"lat"};
573                next if ( ($x2==$x1) && ($y2==$y1) ); # in an ideal world this would be unnecessary, but as we do numerical calculations...
574                next if ( ($x3==$x1) && ($y3==$y1) ); # in an ideal world this would be unnecessary, but as we do numerical calculations...
575                my ($ret)= intersect($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3);
576                if ($ret>0) # found intersection
577                {
578                        printf("   SEG %3d: %d (%f,%f) --> %d (%f,%f) \n", $seg, $from, $x2, $y2, $to, $x3, $y3 ) if ($debug);
579                        return 1;
580                }
581        }
582        return 0; # no intersection
583}
584
585
586# Perl implementation of intersect algorithm for segments
587# first segments is from ($x0,$y0) to ($x1,$y1), second segment is from ($x2,$y2) to ($x3,$y3)
588# source: http://softsurfer.com/Archive/algorithm_0104/algorithm_0104B.htm#Line%20Intersections
589# Modification: This implementation accepts only _true_ intersections, i.e. if the two segments
590# only "touch" at the ends it's _not_ detected as an intersection _intentionally_
591sub intersect
592{
593        my ($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3)= @_;
594       
595        my $cpproduct= ($x1-$x0)*($y3-$y2) - ($y1-$y0)*($x3-$x2);
596        if (abs($cpproduct)>0.000000000001)
597        {
598                my $t1= (($x1-$x0)*($y0-$y2) - ($y1-$y0)*($x0-$x2)) / $cpproduct;
599                return 0 if ( ($t1<=0) || ($t1>=1) );
600                my $t2= (($x3-$x2)*($y0-$y2) - ($y3-$y2)*($x0-$x2)) / $cpproduct;
601                return 0 if ( ($t2<=0) || ($t2>=1) );
602                # a true intersection
603                my $xs= $x0 + $t2*($x1-$x0);
604                my $ys= $y0 + $t2*($y1-$y0);
605                return (1, $xs, $ys);
606        }
607        else
608        {
609                # are the two segments collinear ?
610                my $cpp1= ($x1-$x0)*($y0-$y2) -($y1-$y0)*($x0-$x2);
611                my $cpp2= ($x3-$x2)*($y0-$y2) -($y3-$y2)*($x0-$x2);
612                return 0 if ( ($cpp1!=0) || ($cpp2!=0) ); # they are _not_ collinear
613               
614                return -1 if ( ($x0==$x1) && ($y0==$y1) ); # error: first segment is only a single point
615                return -1 if ( ($x2==$x3) && ($y2==$y3) ); # error: second segment is only a single point
616               
617                my ($t0, $t1);
618                if (($x3-$x2)!=0)
619                {
620                        $t0= ($x0-$x2) / ($x3-$x2);
621                        $t1= ($x1-$x2) / ($x3-$x2);
622                }
623                else
624                {
625                        $t0= ($y0-$y2) / ($y3-$y2);
626                        $t1= ($y1-$y2) / ($y3-$y2);
627                }
628                if ($t0>$t1)
629                {
630                        my $tmp= $t0; $t0= $t1; $t1= $tmp; 
631                }
632                return 0 if ( ($t0>1) || ($t1<0) );
633                $t0 = ($t0<0)? 0 : $t0; # clip to min 0
634                $t1 = ($t1>1)? 1 : $t1; # clip to max 1         
635                if ($t0==$t1)
636                {
637                        return (1, $x2 + $t0*($x3-$x2), $y2 + $t0*($y3-$y2) );
638                }
639                return (2, $x2 + $t0*($x3-$x2), $y2 + $t0*($y3-$y2), $x2 + $t1*($x3-$x2), $y2 + $t1*($y3-$y2));
640        }
641}
642
643
644
645sub make_node
646{
647    my ($lat, $lon) = @_;
648    my $id = --$nodecount;
649    print "<node id='$id' lat='$lat' lon='$lon' />\n";
650    # save node in array for later processing
651    $nodes->{$id}={"lat"=>$lat, "lon"=>$lon};
652    return $id;
653}
654
655sub make_seg
656{
657    my ($from, $to) = @_;
658    my $id = ++$segcount;
659    $segments->{$id} = { "from" => $from, "to" => $to };
660    return $id;
661}
662
663sub make_way
664{
665    my ($seglist, $open) = @_;
666    my $id = --$waycount;
667    print "<way id='$id'>\n";
668    print "  <tag k='natural' v='coastline' />\n";
669    print "  <tag k='created-with' v='close-areas.pl' />\n";
670    print "  <tag k='close-areas.pl:debug' v='open way' />\n" if ($open);
671
672    my $first = 1;
673    foreach my $seg(@$seglist) 
674    { 
675        printf "  <nd ref='%s' />\n", $segments->{$seg}->{"from"} if ($first);
676        printf "  <nd ref='%s' />\n", $segments->{$seg}->{"to"};
677        $first = 0;
678    }
679
680    print "</way>\n";
681    return $id;
682}
683
684sub make_multipolygon
685{
686    my ($seglistlist) = @_;
687    return unless scalar(@$seglistlist);
688
689    # we will create a list of multipolygons. each multipolygon is a list of
690    # ways, with the first being the "outer" polygon.
691    my $multipolygons = [];
692       
693    # first create all the ways.
694    my $way_ids = [];
695    for (my $i=0; $i<scalar(@$seglistlist); $i++)
696    {
697        $way_ids->[$i] = make_way($seglistlist->[$i]);
698    }
699
700    # now find out which way is contained in which. since our polygons do not
701    # intersect, we simply take one node of a way and check whether this node
702    # is contained in any of the other ways; if yes, the whole way is contained.
703    my $contained_in = [];
704    my $children = [];
705    for (my $i=0; $i<scalar(@$seglistlist); $i++)
706    {
707        my $pt = $nodes->{$segments->{$seglistlist->[$i]->[0]}->{"from"}};
708        for (my $j = 0; $j < scalar(@$seglistlist); $j++) # fixme really need to check all?
709        {
710            next if ($j==$i);
711            if (polygon_contains_point($seglistlist->[$j], $pt))
712            {
713                $contained_in->[$i] = $j; 
714                $children->[$j]++;
715                last;
716            }
717        }
718    }
719
720    # now write multipolygons
721    for (my $i=0; $i<scalar(@$seglistlist); $i++)
722    {
723        # ways that don't have children do not trigger a relation
724        next unless $children->[$i];
725        my $id = --$relcount;
726        print "<relation id='$id'>\n";
727        print "  <tag k='type' v='multipolygon' />\n";
728        printf "  <member type='way' ref='%s' role='outer' />\n", $way_ids->[$i];
729        for (my $j = 0; $j < scalar(@$seglistlist); $j++) 
730        {
731            if (defined($contained_in->[$j]) && $contained_in->[$j] == $i)
732            {
733                printf "  <member type='way' ref='%s' role='inner' />\n", $way_ids->[$j];
734            }
735        }
736        print "</relation>\n";
737    }
738}
739
740# index lookup by Martijn van Oosterhout
741sub lookup_handler
742{
743    my ($helpernodes, $x, $y, $zoom) = @_;
744    # make it use z12 x,y coordinates
745    # this looks up the most upper left z12 tile in zoom<12. This probably
746    # needs to be made smarter for islands etc.
747    my $tilex = $x*(2**(12-$zoom));
748    my $tiley = $y*(2**(12-$zoom));
749    my $tileoffset = ($tiley * (2**12)) + $tilex;
750    my $fh;
751    open($fh, "<", "png2tileinfo/oceantiles_12.dat") or die;
752    seek $fh, int($tileoffset / 4), 0; 
753    my $buffer;
754    read $fh, $buffer, 1;
755    $buffer = substr( $buffer."\0", 0, 1 );
756    $buffer = unpack "B*", $buffer;
757    my $str = substr( $buffer, 2*($tileoffset % 4), 2 );
758    close($fh);
759
760    print("lookup handler finds: $str\n") if $debug;
761
762    # $str eq "00" => unknown (not yet checked)
763    # $str eq "01" => known land
764    # $str eq "10" => known sea
765    # $str eq "11" => known edge tile
766
767    return $str; 
768}
769
770sub addBlueRectangle
771{
772    my $helpernodes = shift;
773    my @n;
774    my @s;
775    for (my $i=0; $i<4; $i++)
776    {
777        $n[$i] = make_node($helpernodes->{$i}->[0], 
778                         $helpernodes->{$i}->[1]);
779    }
780    for (my $i=0; $i<4; $i++)
781    {
782        if ($direction eq "ccw")
783        {
784            $s[$i] = make_seg($n[$i], $n[($i+1)%4]);
785        }
786        else
787        {
788            $s[3-$i] = make_seg($n[($i+3)%4], $n[($i+2)%4]);
789        }
790    }
791    push(@coastline_segments, \@s);
792}
793
794# this takes two points (hashes with lat/lon keys) as arguments and returns
795# an array reference to an array containing up to four points (hashes with
796# lat/lon keys and an added "side" key, where 0=right 1=top 2=left 3=bottom)
797# denoting the intersections of the line formed by the input points with the
798# bounding box.
799#
800# 0 results - segment is fully inside or fully outside bounding box
801# 1 result  - segment begins inside, ends outside or vice versa
802# 2 results - segment begins outside and ends outside, but cuts through
803#             bounding box
804# 3 results - can't think how this can happen
805# 4 results - as case 2 but segment is exactly at the bounding box diagonal,
806#             thus intersecting with all four edges (freak case)
807#
808# If there is more than one result, results are sorted by distance from the
809# first point.
810sub compute_bbox_intersections
811{
812    my ($from, $to) = @_;
813    my $result = [];
814    my $point;
815
816    # div by zero FIXME
817    my $latd = ($to->{"lat"} - $from->{"lat"});
818    my $lond = ($to->{"lon"} - $from->{"lon"});
819
820    # segment's bbox
821    my $s_minlon = $from->{"lon"};
822    my $s_minlat = $from->{"lat"};
823    $s_minlon = $to->{"lon"} if ($to->{"lon"} < $s_minlon);
824    $s_minlat = $to->{"lat"} if ($to->{"lat"} < $s_minlat);
825    my $s_maxlon = $from->{"lon"};
826    my $s_maxlat = $from->{"lat"};
827    $s_maxlon = $to->{"lon"} if ($to->{"lon"} > $s_maxlon);
828    $s_maxlat = $to->{"lat"} if ($to->{"lat"} > $s_maxlat);
829
830    printf "BBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n",
831        $minlat, $minlon, $maxlat, $maxlon if ($debug);
832
833    printf "SBBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n",
834        $s_minlat, $s_minlon, $s_maxlat, $s_maxlon if ($debug);
835
836    # only if the segment is not horizontal
837    if ($latd != 0)
838    {
839        # intersection with top of bounding box
840        $point = { 
841            "side" => "1",
842            "lat" => $maxlat, 
843            "lon" => $from->{"lon"} + ($maxlat - $from->{"lat"}) * $lond / $latd 
844        };
845        push (@$result, $point) 
846            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
847               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
848               ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon);
849
850
851        # intersection with bottom of bounding box
852        $point = { 
853            "side" => "3",
854            "lat" => $minlat, 
855            "lon" => $from->{"lon"} + ($minlat - $from->{"lat"}) * $lond / $latd 
856        };
857        push (@$result, $point) 
858            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
859               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
860               ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon);
861    }
862
863    # only if the segment is not vertical
864    if ($lond != 0)
865    {
866        # intersection with left of bounding box
867        $point = { 
868            "side" => "2",
869            "lat" => $from->{"lat"} + $latd / $lond * ($minlon - $from->{"lon"}),
870            "lon" => $minlon
871        };
872        push (@$result, $point) 
873            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
874               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
875               ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat);
876
877        # intersection with right of bounding box
878        $point = { 
879            "side" => "0",
880            "lat" => $from->{"lat"} + $latd / $lond * ($maxlon - $from->{"lon"}),
881            "lon" => $maxlon
882        };
883        push (@$result, $point) 
884            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
885               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
886               ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat);
887    }
888
889    # if more than 1 result, sort by distance from origin of segment
890    # (strictly speaking this sorts by distance squared but why waste
891    # a sqrt call)
892    if (scalar(@$result) > 1)
893    {
894        my @tmp = sort 
895            { (($a->{"lat"} - $from->{"lat"})**2 + ($a->{"lon"} - $from->{"lon"})**2) <=> 
896              (($b->{"lat"} - $from->{"lat"})**2 + ($b->{"lon"} - $from->{"lon"})**2) } @$result;
897        $result = \@tmp;
898    }
899
900    #printf "intersections for segment %f,%f - %f,%f:\n", $from->{"lat"}, $from->{"lon"}, $to->{"lat"}, $to->{"lon"};
901    #print Dumper($result);
902
903    return $result;
904}
905
906# expects point as lat/lon hash, and polygon as a list of segment ids
907# which will be looked up in global $segments to get from/to node ids
908# which will in turn be looked up in $nodes to get lat/lon
909sub polygon_contains_point
910{
911    my ($seglist, $point) = @_;
912    my $p1 =  $nodes->{$segments->{$seglist->[0]}->{"from"}};
913    my $counter = 0;
914    foreach my $seg(@$seglist)
915    {
916        my $p2 = $nodes->{$segments->{$seg}->{"to"}};
917        if ($point->{"lat"} >= $p1->{"lat"} || $point->{"lat"} >= $p2->{"lat"})
918        {
919            if ($point->{"lat"} < $p1->{"lat"} || $point->{"lat"} < $p2->{"lat"})
920            {
921                if ($point->{"lon"} < $p1->{"lon"} || $point->{"lon"} < $p2->{"lon"})
922                {
923                    if ($p1->{"lat"} != $p2->{"lat"})
924                    {
925                        my $xint = ($point->{"lat"}-$p1->{"lat"})*($p2->{"lon"}-$p1->{"lon"})/($p2->{"lat"}-$p1->{"lat"})+$p1->{"lon"};
926                        if ($p1->{"lon"} == $p2->{"lon"} || $point->{"lon"} <= $xint)
927                        {
928                            $counter++;
929                        }
930                    }
931                }
932            }
933        }
934        $p1 = $p2;
935    }
936    return ($counter%2);
937}
938
939sub node_is_inside
940{
941    my $point = shift;
942    return 0 if ($point->{"lat"} > $maxlat);
943    return 0 if ($point->{"lat"} < $minlat);
944    return 0 if ($point->{"lon"} > $maxlon);
945    return 0 if ($point->{"lon"} < $minlon);
946    return 1;
947}
948
949
950sub compute_angle_from_bbox_center
951{
952    my ($node) = @_;
953    my $opposite_leg = ($node->{"lat"}-($maxlat-$minlat)/2-$minlat);
954    my $adjacent_leg = ($node->{"lon"}-($maxlon-$minlon)/2-$minlon);
955    my $z = cplx($adjacent_leg, $opposite_leg);
956    return (arg($z) < 0) ? arg($z) + 2*pi : arg($z);
957}
958
Note: See TracBrowser for help on using the repository browser.