source: subversion/applications/rendering/tilesAtHome/close-areas.pl @ 15164

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

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

File size: 30.7 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");
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            my $newnode = make_node($helpernodes->{$i}->[0], $helpernodes->{$i}->[1]);
447            my $newseg = make_seg($currentnode, $newnode);
448            $currentnode = $newnode;
449            push(@seglist, $newseg);
450        }
451    }
452    else
453    {
454        my $min_once;
455        if ($exit_side == 0 and $entry_side == 0)
456        {
457            $min_once = (($exit_angle < $entry_angle and ($entry_angle - $exit_angle) < pi()) ||
458                         ($exit_angle > pi() and $entry_angle < pi()));
459        }
460        else
461        {
462            $min_once = ($exit_angle < $entry_angle);
463        }
464        printf("min_once=%d\n", $min_once) if $debug;
465        for (my $i = $exit_side;; $i--)
466        {
467            $i=3 if ($i==-1);
468            last if ($i==$entry_side && !$min_once);
469            $min_once = 0;
470            my $helper = $i-1;
471            $helper = 3 if ($helper == -1);
472            last if (check_segments_intersect($helpernodes->{$helper}->[1], $helpernodes->{$helper}->[0], $nodes->{$currentnode}->{"lon"}, $nodes->{$currentnode}->{"lat"}, $segments) ); # new
473            my $newnode = make_node($helpernodes->{$helper}->[0], $helpernodes->{$helper}->[1]);
474            my $newseg = make_seg($currentnode, $newnode);
475            $currentnode = $newnode;
476            push(@seglist, $newseg);
477        }
478    }
479
480    my $newseg = make_seg($currentnode, $segments->{$found}->{"from"});
481    push(@seglist, $newseg);
482
483    # if the segment we have found is our start segment (which we added
484    # as a special case earlier!), we have a closed way.
485    if ($found == $seglist[0])
486    {
487        printf("CLOSED\n") if ($debug);
488        push(@coastline_segments, \@seglist);
489        next;
490    }
491
492    # else we just plod on.
493    push (@seglist, $found);
494    $currentseg = $found;
495    $currentnode = $segments->{$found}->{"to"};
496    delete $available_segs->{$found};
497    #nprintf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
498    goto STRING;
499}
500
501ENDOSM:
502
503# if we had no coastline intersection with the bounding box, but coastline
504# was present, we have an island situation and need to add a blue background.
505# FIXME what if we have a "little lake" situation?
506
507unless( $border_crossed )
508{
509  my $state = lookup_handler($helpernodes, $tilex, $tiley, $zoom);
510  if( $state eq "10" )
511  {
512    # sea
513    addBlueRectangle($helpernodes);
514  }
515  elsif ( $state eq "01" )
516  {
517    #land
518  }
519  else # state 00 or 11 (unknown or mixed)
520  { 
521    my %temp = ("00"=>0, "10"=>0, "01"=>0, "11"=>0);;
522    $temp{lookup_handler($helpernodes, $tilex-1, $tiley, $zoom)}++;
523    $temp{lookup_handler($helpernodes, $tilex+1, $tiley, $zoom)}++;
524    $temp{lookup_handler($helpernodes, $tilex, $tiley-1, $zoom)}++;
525    $temp{lookup_handler($helpernodes, $tilex, $tiley+1, $zoom)}++;
526
527    if( $temp{"10"} + $temp{"11"} > $temp{"01"} ) # if more sea/coast than land around.
528    {
529      addBlueRectangle($helpernodes);
530    }
531    elsif ( ($state eq "11") and ( $temp{"01"} == 0 ) ) 
532    # if the tile is marked coast but no land near, assume it's a group of islands instead of lakes.
533    {
534      # coast
535      addBlueRectangle($helpernodes);
536    }
537    else
538    {
539      #land
540    }
541  }
542}
543
544if (scalar(@coastline_segments) == 1)
545{
546    make_way($coastline_segments[0]);
547}
548else
549{
550    make_multipolygon(\@coastline_segments);
551}
552print "</osm>\n";
553
554
555
556
557# checks segment from ($x0,$y0) to ($x1,$y1) for intersection with all segments
558# in $s
559sub check_segments_intersect
560{
561        my ($x0, $y0, $x1, $y1, $s)= @_;
562       
563        printf("check_segments_intersect: (%f,%f)  (%f,%f)   %s\n", $x0, $y0, $x1, $y1, join(", ",keys(%$s)) ) if ($debug);
564        foreach my $seg(keys(%$s))
565        {
566                my $from= $s->{$seg}->{"from"};
567                my $to  = $s->{$seg}->{"to"};
568                my $x2= $nodes->{$from}->{"lon"};
569                my $y2= $nodes->{$from}->{"lat"};
570                my $x3= $nodes->{$to}  ->{"lon"};
571                my $y3= $nodes->{$to}  ->{"lat"};
572                my ($ret)= intersect($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3);
573                if ($ret>0) # found intersection
574                {
575                        printf("   SEG %3d: %d (%f,%f) --> %d (%f,%f) \n", $seg, $from, $x2, $y2, $to, $x3, $y3 ) if ($debug);
576                        return 1;
577                }
578        }
579        return 0; # no intersection
580}
581
582
583# Perl implementation of intersect algorithm for segments
584# first segments is from ($x0,$y0) to ($x1,$y1), second segment is from ($x2,$y2) to ($x3,$y3)
585# source: http://softsurfer.com/Archive/algorithm_0104/algorithm_0104B.htm#Line%20Intersections
586# Modification: This implementation accepts only _true_ intersections, i.e. if the two segments
587# only "touch" at the ends it's _not_ detected as an intersection _intentionally_
588sub intersect
589{
590        my ($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3)= @_;
591       
592        my $cpproduct= ($x1-$x0)*($y3-$y2) - ($y1-$y0)*($x3-$x2);
593        if (abs($cpproduct)>0.000000000001)
594        {
595                my $t1= (($x1-$x0)*($y0-$y2) - ($y1-$y0)*($x0-$x2)) / $cpproduct;
596                return 0 if ( ($t1<=0) || ($t1>=1) );
597                my $t2= (($x3-$x2)*($y0-$y2) - ($y3-$y2)*($x0-$x2)) / $cpproduct;
598                return 0 if ( ($t2<=0) || ($t2>=1) );
599                # a true intersection
600                my $xs= $x0 + $t2*($x1-$x0);
601                my $ys= $y0 + $t2*($y1-$y0);
602                return (1, $xs, $ys);
603        }
604        else
605        {
606                # are the two segments collinear ?
607                my $cpp1= ($x1-$x0)*($y0-$y2) -($y1-$y0)*($x0-$x2);
608                my $cpp2= ($x3-$x2)*($y0-$y2) -($y3-$y2)*($x0-$x2);
609                return 0 if ( ($cpp1!=0) || ($cpp2!=0) ); # they are _not_ collinear
610               
611                return -1 if ( ($x0==$x1) && ($y0==$y1) ); # error: first segment is only a single point
612                return -1 if ( ($x2==$x3) && ($y2==$y3) ); # error: second segment is only a single point
613               
614                my ($t0, $t1);
615                if (($x3-$x2)!=0)
616                {
617                        $t0= ($x0-$x2) / ($x3-$x2);
618                        $t1= ($x1-$x2) / ($x3-$x2);
619                }
620                else
621                {
622                        $t0= ($y0-$y2) / ($y3-$y2);
623                        $t1= ($y1-$y2) / ($y3-$y2);
624                }
625                if ($t0>$t1)
626                {
627                        my $tmp= $t0; $t0= $t1; $t1= $tmp; 
628                }
629                return 0 if ( ($t0>1) || ($t1<0) );
630                $t0 = ($t0<0)? 0 : $t0; # clip to min 0
631                $t1 = ($t1>1)? 1 : $t1; # clip to max 1         
632                if ($t0==$t1)
633                {
634                        return (1, $x2 + $t0*($x3-$x2), $y2 + $t0*($y3-$y2) );
635                }
636                return (2, $x2 + $t0*($x3-$x2), $y2 + $t0*($y3-$y2), $x2 + $t1*($x3-$x2), $y2 + $t1*($y3-$y2));
637        }
638}
639
640
641
642sub make_node
643{
644    my ($lat, $lon) = @_;
645    my $id = --$nodecount;
646    print "<node id='$id' lat='$lat' lon='$lon' />\n";
647    # save node in array for later processing
648    $nodes->{$id}={"lat"=>$lat, "lon"=>$lon};
649    return $id;
650}
651
652sub make_seg
653{
654    my ($from, $to) = @_;
655    my $id = ++$segcount;
656    $segments->{$id} = { "from" => $from, "to" => $to };
657    return $id;
658}
659
660sub make_way
661{
662    my ($seglist, $open) = @_;
663    my $id = --$waycount;
664    print "<way id='$id'>\n";
665    print "  <tag k='natural' v='coastline' />\n";
666    print "  <tag k='created-with' v='close-areas.pl' />\n";
667    print "  <tag k='close-areas.pl:debug' v='open way' />\n" if ($open);
668
669    my $first = 1;
670    foreach my $seg(@$seglist) 
671    { 
672        printf "  <nd ref='%s' />\n", $segments->{$seg}->{"from"} if ($first);
673        printf "  <nd ref='%s' />\n", $segments->{$seg}->{"to"};
674        $first = 0;
675    }
676
677    print "</way>\n";
678    return $id;
679}
680
681sub make_multipolygon
682{
683    my ($seglistlist) = @_;
684    return unless scalar(@$seglistlist);
685
686    # we will create a list of multipolygons. each multipolygon is a list of
687    # ways, with the first being the "outer" polygon.
688    my $multipolygons = [];
689       
690    # first create all the ways.
691    my $way_ids = [];
692    for (my $i=0; $i<scalar(@$seglistlist); $i++)
693    {
694        $way_ids->[$i] = make_way($seglistlist->[$i]);
695    }
696
697    # now find out which way is contained in which. since our polygons do not
698    # intersect, we simply take one node of a way and check whether this node
699    # is contained in any of the other ways; if yes, the whole way is contained.
700    my $contained_in = [];
701    my $children = [];
702    for (my $i=0; $i<scalar(@$seglistlist); $i++)
703    {
704        my $pt = $nodes->{$segments->{$seglistlist->[$i]->[0]}->{"from"}};
705        for (my $j = 0; $j < scalar(@$seglistlist); $j++) # fixme really need to check all?
706        {
707            next if ($j==$i);
708            if (polygon_contains_point($seglistlist->[$j], $pt))
709            {
710                $contained_in->[$i] = $j; 
711                $children->[$j]++;
712                last;
713            }
714        }
715    }
716
717    # now write multipolygons
718    for (my $i=0; $i<scalar(@$seglistlist); $i++)
719    {
720        # ways that don't have children do not trigger a relation
721        next unless $children->[$i];
722        my $id = --$relcount;
723        print "<relation id='$id'>\n";
724        print "  <tag k='type' v='multipolygon' />\n";
725        printf "  <member type='way' ref='%s' role='outer' />\n", $way_ids->[$i];
726        for (my $j = 0; $j < scalar(@$seglistlist); $j++) 
727        {
728            if (defined($contained_in->[$j]) && $contained_in->[$j] == $i)
729            {
730                printf "  <member type='way' ref='%s' role='inner' />\n", $way_ids->[$j];
731            }
732        }
733        print "</relation>\n";
734    }
735}
736
737# index lookup by Martijn van Oosterhout
738sub lookup_handler
739{
740    my ($helpernodes, $x, $y, $zoom) = @_;
741    # make it use z12 x,y coordinates
742    # this looks up the most upper left z12 tile in zoom<12. This probably
743    # needs to be made smarter for islands etc.
744    my $tilex = $x*(2**(12-$zoom));
745    my $tiley = $y*(2**(12-$zoom));
746    my $tileoffset = ($tiley * (2**12)) + $tilex;
747    my $fh;
748    open($fh, "<", "png2tileinfo/oceantiles_12.dat") or die;
749    seek $fh, int($tileoffset / 4), 0; 
750    my $buffer;
751    read $fh, $buffer, 1;
752    $buffer = substr( $buffer."\0", 0, 1 );
753    $buffer = unpack "B*", $buffer;
754    my $str = substr( $buffer, 2*($tileoffset % 4), 2 );
755    close($fh);
756
757    print("lookup handler finds: $str\n") if $debug;
758
759    # $str eq "00" => unknown (not yet checked)
760    # $str eq "01" => known land
761    # $str eq "10" => known sea
762    # $str eq "11" => known edge tile
763
764    return $str; 
765}
766
767sub addBlueRectangle
768{
769    my $helpernodes = shift;
770    my @n;
771    my @s;
772    for (my $i=0; $i<4; $i++)
773    {
774        $n[$i] = make_node($helpernodes->{$i}->[0], 
775                         $helpernodes->{$i}->[1]);
776    }
777    for (my $i=0; $i<4; $i++)
778    {
779        if ($direction eq "ccw")
780        {
781            $s[$i] = make_seg($n[$i], $n[($i+1)%4]);
782        }
783        else
784        {
785            $s[3-$i] = make_seg($n[($i+3)%4], $n[($i+2)%4]);
786        }
787    }
788    push(@coastline_segments, \@s);
789}
790
791# this takes two points (hashes with lat/lon keys) as arguments and returns
792# an array reference to an array containing up to four points (hashes with
793# lat/lon keys and an added "side" key, where 0=right 1=top 2=left 3=bottom)
794# denoting the intersections of the line formed by the input points with the
795# bounding box.
796#
797# 0 results - segment is fully inside or fully outside bounding box
798# 1 result  - segment begins inside, ends outside or vice versa
799# 2 results - segment begins outside and ends outside, but cuts through
800#             bounding box
801# 3 results - can't think how this can happen
802# 4 results - as case 2 but segment is exactly at the bounding box diagonal,
803#             thus intersecting with all four edges (freak case)
804#
805# If there is more than one result, results are sorted by distance from the
806# first point.
807sub compute_bbox_intersections
808{
809    my ($from, $to) = @_;
810    my $result = [];
811    my $point;
812
813    # div by zero FIXME
814    my $latd = ($to->{"lat"} - $from->{"lat"});
815    my $lond = ($to->{"lon"} - $from->{"lon"});
816
817    # segment's bbox
818    my $s_minlon = $from->{"lon"};
819    my $s_minlat = $from->{"lat"};
820    $s_minlon = $to->{"lon"} if ($to->{"lon"} < $s_minlon);
821    $s_minlat = $to->{"lat"} if ($to->{"lat"} < $s_minlat);
822    my $s_maxlon = $from->{"lon"};
823    my $s_maxlat = $from->{"lat"};
824    $s_maxlon = $to->{"lon"} if ($to->{"lon"} > $s_maxlon);
825    $s_maxlat = $to->{"lat"} if ($to->{"lat"} > $s_maxlat);
826
827    printf "BBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n",
828        $minlat, $minlon, $maxlat, $maxlon if ($debug);
829
830    printf "SBBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n",
831        $s_minlat, $s_minlon, $s_maxlat, $s_maxlon if ($debug);
832
833    # only if the segment is not horizontal
834    if ($latd != 0)
835    {
836        # intersection with top of bounding box
837        $point = { 
838            "side" => "1",
839            "lat" => $maxlat, 
840            "lon" => $from->{"lon"} + ($maxlat - $from->{"lat"}) * $lond / $latd 
841        };
842        push (@$result, $point) 
843            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
844               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
845               ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon);
846
847
848        # intersection with bottom of bounding box
849        $point = { 
850            "side" => "3",
851            "lat" => $minlat, 
852            "lon" => $from->{"lon"} + ($minlat - $from->{"lat"}) * $lond / $latd 
853        };
854        push (@$result, $point) 
855            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
856               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
857               ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon);
858    }
859
860    # only if the segment is not vertical
861    if ($lond != 0)
862    {
863        # intersection with left of bounding box
864        $point = { 
865            "side" => "2",
866            "lat" => $from->{"lat"} + $latd / $lond * ($minlon - $from->{"lon"}),
867            "lon" => $minlon
868        };
869        push (@$result, $point) 
870            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
871               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
872               ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat);
873
874        # intersection with right of bounding box
875        $point = { 
876            "side" => "0",
877            "lat" => $from->{"lat"} + $latd / $lond * ($maxlon - $from->{"lon"}),
878            "lon" => $maxlon
879        };
880        push (@$result, $point) 
881            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
882               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
883               ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat);
884    }
885
886    # if more than 1 result, sort by distance from origin of segment
887    # (strictly speaking this sorts by distance squared but why waste
888    # a sqrt call)
889    if (scalar(@$result) > 1)
890    {
891        my @tmp = sort 
892            { (($a->{"lat"} - $from->{"lat"})**2 + ($a->{"lon"} - $from->{"lon"})**2) <=> 
893              (($b->{"lat"} - $from->{"lat"})**2 + ($b->{"lon"} - $from->{"lon"})**2) } @$result;
894        $result = \@tmp;
895    }
896
897    #printf "intersections for segment %f,%f - %f,%f:\n", $from->{"lat"}, $from->{"lon"}, $to->{"lat"}, $to->{"lon"};
898    #print Dumper($result);
899
900    return $result;
901}
902
903# expects point as lat/lon hash, and polygon as a list of segment ids
904# which will be looked up in global $segments to get from/to node ids
905# which will in turn be looked up in $nodes to get lat/lon
906sub polygon_contains_point
907{
908    my ($seglist, $point) = @_;
909    my $p1 =  $nodes->{$segments->{$seglist->[0]}->{"from"}};
910    my $counter = 0;
911    foreach my $seg(@$seglist)
912    {
913        my $p2 = $nodes->{$segments->{$seg}->{"to"}};
914        if ($point->{"lat"} >= $p1->{"lat"} || $point->{"lat"} >= $p2->{"lat"})
915        {
916            if ($point->{"lat"} < $p1->{"lat"} || $point->{"lat"} < $p2->{"lat"})
917            {
918                if ($point->{"lon"} < $p1->{"lon"} || $point->{"lon"} < $p2->{"lon"})
919                {
920                    if ($p1->{"lat"} != $p2->{"lat"})
921                    {
922                        my $xint = ($point->{"lat"}-$p1->{"lat"})*($p2->{"lon"}-$p1->{"lon"})/($p2->{"lat"}-$p1->{"lat"})+$p1->{"lon"};
923                        if ($p1->{"lon"} == $p2->{"lon"} || $point->{"lon"} <= $xint)
924                        {
925                            $counter++;
926                        }
927                    }
928                }
929            }
930        }
931        $p1 = $p2;
932    }
933    return ($counter%2);
934}
935
936sub node_is_inside
937{
938    my $point = shift;
939    return 0 if ($point->{"lat"} > $maxlat);
940    return 0 if ($point->{"lat"} < $minlat);
941    return 0 if ($point->{"lon"} > $maxlon);
942    return 0 if ($point->{"lon"} < $minlon);
943    return 1;
944}
945
946
947sub compute_angle_from_bbox_center
948{
949    my ($node) = @_;
950    my $opposite_leg = ($node->{"lat"}-($maxlat-$minlat)/2-$minlat);
951    my $adjacent_leg = ($node->{"lon"}-($maxlon-$minlon)/2-$minlon);
952    my $z = cplx($adjacent_leg, $opposite_leg);
953    return (arg($z) < 0) ? arg($z) + 2*pi : arg($z);
954}
955
Note: See TracBrowser for help on using the repository browser.