source: subversion/applications/utils/osm-extract/polygons/extract-polygon.pl

Last change on this file was 5345, checked in by frederik, 12 years ago

0.5 versions now up to date

  • Property svn:executable set to *
File size: 18.2 KB
Line 
1#!/usr/bin/perl
2
3# This script extracts all nodes, segments, and ways lying inside a polygon
4#
5# The script first compiles a polygon from the data, then simplifies it,
6# and then processes the planet file. In processing the file, first a
7# bounding box is used to exclude all nodes that lie outside, and the
8# more expensive polygon check is only made for those inside.
9#
10# Author Frederik Ramm <frederik@remote.org> / public domain
11# Author Keith Sharp <kms@passback.co.uk> / public domain
12
13# Progress meter, exclusions, et al all added by Martijn van Oosterhout
14# <kleptog@svana.org> for the AND import (August 2007) / public domain
15
16# THIS IS THE API 0.5 VERSION (that's 0.4 minus segments plus relations).
17# Also adapted to used Bit::Vector instead of hash to save memory. Big
18# change is that this now supports "referential integrity", i.e. you
19# can request it to behave like the API does, giving you all nodes for
20# each way returned, not only those within the bounding box.
21#
22# Get this from CPAN!
23use Math::Polygon;
24use Getopt::Long;
25use File::Basename;
26use Bit::Vector;
27use IO::Handle;
28use IO::File;
29
30use strict;
31
32my $borderpolys;
33my $currentpoints;
34my $minlon = 999;
35my $minlat = 999;
36my $maxlat = -999;
37my $maxlon = -999;
38
39# assume 500 million nodes,
40my $nodes_copied = Bit::Vector->new(500*1000*1000);
41my $nodes_needed = Bit::Vector->new(500*1000*1000);
42
43# and 50 million ways and relations.
44my $ways_copied = Bit::Vector->new(50*1000*1000);
45my $ways_needed = Bit::Vector->new(50*1000*1000);
46
47my $relations_copied = Bit::Vector->new(50*1000*1000);
48my $relations_needed = Bit::Vector->new(50*1000*1000);
49
50# this puts total memory allocation at 1,2 billion bits, which should
51# be just below 250 MB including overhead and should be sufficient for
52# planet files in the area of 5 to 10 GB (compressed size).
53
54my $help;
55my $polyfile;
56my $infile;
57my $outfile;
58my $remainsfile;
59my $compress;
60my $verbose = 0;
61my $references = 0;
62
63my $prog = basename($0);
64
65my $filepos_nodes = 1;
66my $filepos_ways;
67my $filepos_relations;
68
69my $input_file_handle;
70my $input_is_seekable;
71my $assume_input_sorted = 1;
72my $deep_relations = 1;
73my $sort_output = 0;
74my $bbox;
75my $area = 0;
76
77my $output_file_handle_nodes;
78my $output_file_handle_ways;
79my $output_file_handle_rels;
80
81sub usage 
82{
83    my ($msg) = @_;
84    print STDERR "$msg\n" if defined($msg);
85
86    print STDERR << "EOF";
87
88This Perl script will process a planet.osm file and extract the nodes,
89ways, and relations falling within a polygon.
90
91usage: $prog [-h] [-a number] [-c number] [-i file] [-o file] [-r file] -p file [-d]
92
93 -h       : print ths help message and exit.
94 -a num   : ignore sub-polygons of area "num" or smaller (in degrees squared)
95 -c num   : reduce the number of points in the polygon by num %.  This 
96            only takes effect if the polgon has more than 100 points.
97 -b bbox  : use bounding box instead of polygon (bbox being four comma-
98            separated numbers: bllon,bllat,trlon,trlat).
99 -p file  : file containing the polygon definition.
100 -i file  : OSM planet file to process.
101 -o file  : OSM planet file to output.
102 -r       : preserve referential integrity. This requires multiple passes
103            and thus will be slower, ESPECIALLY if you use a compressed
104            input file (needs to be decompressed multiple times). 
105 -s       : sort output; this is only meaningful in conjunction with -r
106            as output without -r will be sorted anyway.
107
108If you do not supply an output file, will write to STDOUT. Reading from
109STDIN is not yet supported.
110
111A polygon file should be a plain text file structured like this: Each
112polygon begins with a number in the first column; the following rows
113begin with whitespace followed by a whitespace-separated coordinate
114pair (longitude, then latitude). The polygon ends when a line contains
115only "END". Other polygons may then follow.
116
117For example a polygon defining Great Britain and Ireland (and some of the
118smaller islands) would be:
119
1201
121    -0.6450E+01    0.4980E+02
122    -0.2000E+01    0.4890E+02
123    -0.1850E+01    0.4925E+02
124    -0.2080E+01    0.4973E+02
125     0.1350E+01    0.5090E+02
126     0.2250E+01    0.5258E+02
127    -0.0500E+01    0.6130E+02
128    -0.8920E+01    0.5785E+02
129    -0.1140E+02    0.5130E+02
130    -0.6450E+01    0.4980E+02
131END
132
133Multiple polygons can be described by adding additional sections to the
134polygon file.
135
136EOF
137    exit;
138}
139
140GetOptions ('h|help' => \$help,
141            'c|compress=i' => \$compress,
142            'a|area=n' => \$area,
143            'v|verbose' => \$verbose,
144            'i|input=s' => \$infile,
145            'o|output=s' => \$outfile,
146            'p|polygon=s' => \$polyfile,
147            'b|bbox=s' => \$bbox,
148            's|sort' => \$sort_output,
149            'r|references' => \$references,
150            ) || usage ();
151
152if ($help) {
153    usage();
154}
155
156if ((!$polyfile && !$bbox) || ($polyfile && $bbox))
157{
158    usage("exactly one of -b and -p must be given");
159}
160
161if ($sort_output && !$outfile)
162{
163    usage("you must specify an output file if you want sorting");
164}
165
166if (!$infile)
167{
168    usage("you must specify an input file - cannot work with stdin");
169}
170
171if ($polyfile)
172{
173    $borderpolys = [];
174    open (PF, "$polyfile") || die "Could not open file: $polyfile: $!";
175
176    my $invert;
177    # initialize border polygon.
178    while(<PF>)
179    {
180        if (/^(!?)\d/)
181        {
182            $invert = ($1 eq "!") ? 1 : 0;
183            $currentpoints = [];
184        }
185        elsif (/^END/)
186        {
187            my $pol = Math::Polygon->new(points => $currentpoints);
188            if (($compress > 0 && $compress < 100) && $pol->nrPoints > 99) {
189                my $simp = $pol->simplify($pol->nrPoints*(100-$compress)/100);
190                push(@{$borderpolys}, [$simp,$invert,[$simp->bbox]]) if ($simp->area()>$area);
191            } else {
192            push(@{$borderpolys}, [$pol,$invert,[$pol->bbox]]);
193            }
194            undef $currentpoints;
195            if( $verbose )
196            {
197                printf STDERR "Added polygon: %d points (%.2f,%.2f)-(%.2f,%.2f) %s\n",
198                    $borderpolys->[-1][0]->nrPoints,
199                    @{$borderpolys->[-1][2]},
200                    ($borderpolys->[-1][1] ? "exclude" : "include");
201            }
202        }
203        elsif (defined($currentpoints))
204        {
205            /^\s+([0-9.E+-]+)\s+([0-9.E+-]+)/ or die;
206            push(@{$currentpoints}, [$2, $1]);
207            $minlat = $2 if ($2 < $minlat);
208            $maxlat = $2 if ($2 > $maxlat);
209            $minlon = $1 if ($1 < $minlon);
210            $maxlon = $1 if ($1 > $maxlon);
211        }
212    }
213    close (PF);
214}
215else
216{
217    # simply set minlat, maxlat etc from bbox parameter - no polygons
218    ($minlat, $minlon, $maxlat, $maxlon) = split(",", $bbox);
219    die ("badly formed bounding box - use four comma-separated values for ".
220        "bottom left latitude, bottom left longitude, top right latitude, ".
221        "top right longitude") unless defined($maxlon);
222    die ("max longitude is less than min longitude") if ($maxlon < $minlon);
223    die ("max latitude is less than min latitude") if ($maxlat < $minlat);
224}
225
226if ($outfile) 
227{
228    if ($sort_output)
229    {
230        $output_file_handle_nodes = new IO::File(">$outfile.n") or die ("cannot open $outfile.n: $!");
231        $output_file_handle_ways = new IO::File(">$outfile.w") or die ("cannot open $outfile.w: $!");
232        $output_file_handle_rels = new IO::File(">$outfile.r") or die ("cannot open $outfile.r: $!");
233    }
234    else
235    {
236        $output_file_handle_nodes = new IO::File(">$outfile") or die ("cannot open $outfile: $!");
237        $output_file_handle_ways = $output_file_handle_nodes;
238        $output_file_handle_rels = $output_file_handle_nodes;
239    }
240} 
241else 
242{
243    $output_file_handle_nodes = new IO::Handle;
244    $output_file_handle_nodes->fdopen(fileno(STDOUT),"w") || die "Cannot write to standard output: $!";
245    $output_file_handle_ways = $output_file_handle_nodes;
246    $output_file_handle_rels = $output_file_handle_nodes;
247}
248
249if ($remainsfile) 
250{
251    open (RF, ">$remainsfile") || die "Could not open file: $remainsfile: $!";
252}
253
254print $output_file_handle_nodes <<EOF;
255<?xml version="1.0" encoding="UTF-8"?>
256<osm version="0.5" generator="extract-polygon.pl">
257EOF
258
259my $copy;
260my $waybuf;
261my $count = 0;
262
263my $polygon_node_selector = sub {
264    my ($id, $lat, $lon) = @_;
265    return 0 if (($lat < $minlat) || ($lat > $maxlat));
266    return 0 if (($lon < $minlon) || ($lon > $maxlon));
267    return 1 unless defined($borderpolys);
268    my $ll = [$lat, $lon];
269    my $rv = 0;
270    foreach my $p (@{$borderpolys})
271    {
272        my($poly,$invert,$bbox) = @$p;
273        next if ($ll->[0] < $bbox->[0]) or ($ll->[0] > $bbox->[2]);
274        next if ($ll->[1] < $bbox->[1]) or ($ll->[1] > $bbox->[3]);
275
276        if ($poly->contains($ll))
277        {
278            # If this polygon is for exclusion, we immediately bail and go for the next point
279            if($invert)
280            {
281                return 0;
282            }
283            $rv = 1; 
284            # do not exit here as an exclusion poly may still be there
285        }
286    }
287    return $rv;
288};
289
290# ---------------------------------------------------------------------------
291# First Extraction
292# ---------------------------------------------------------------------------
293
294# This starts copying stuff from the planet file to output, recording IDs
295# as it goes along.
296
297# first, the nodes;
298# this returns a list of nodes having been copied in $used_nodes
299select_nodes($polygon_node_selector, $nodes_copied);
300
301# then, the ways;
302
303my $way_selector = sub {
304    my ($way_id, $node_id) = @_;
305    return ($nodes_copied->contains($node_id));
306};
307
308# this returns a list of ways having been copied in $used_ways,
309# and additionaly a list of all nodes these ways use in $nodes_needed
310
311select_ways($way_selector, $ways_copied, $nodes_needed);
312
313# finally, the relations;
314
315my $relation_selector = sub {
316    my ($relation_id, $member_type, $member_id) = @_;
317    return 1 if ($member_type eq "node" && $nodes_copied->contains($member_id));
318    return 1 if ($member_type eq "way" && $ways_copied->contains($member_id));
319    return 0;
320};
321
322# returns a list of relations having been copied in $relations_copied,
323# and additionally lists of nodes,ways,relations required.
324
325select_relations($relation_selector, $relations_copied);
326
327# use this if you want to allow relations to request the copying of their members:
328# select_relations($relation_selector, $relations_copied, $nodes_needed, $ways_needed, $relations_needed);
329
330# ok, time to pause.
331#
332# what we now have is:
333# * all nodes in the bounding box
334# * all ways using at least one of these nodes
335# * all relations having one of the nodes or ways as their member.
336#
337# That's sufficient for non-referential-integrity mode:
338
339end_output() unless $references;
340
341# ---------------------------------------------------------------------------
342# Follow-On Extractions
343# ---------------------------------------------------------------------------
344
345# if we want referential integrity, we have to re-visit everything, this
346# time collecting referenced objects not collected the first time.
347
348my $node_selector = sub {
349    my ($id, $lat, $lon) = @_;
350    return $nodes_needed->contains($id) && !$nodes_copied->contains($id);
351};
352
353select_nodes($node_selector, $nodes_copied);
354
355# if you want to allow relations to request copying of their members, proceed
356# like this:
357
358#my $way_selector = sub {
359#    my ($way_id, $node_id) = @_;
360#    return 0 if ($ways_copied->contains($way_id));
361#    # this could be used to return all ways that use the copied nodes but
362#    # are not themselves copied:
363#    # return 1 if $nodes_copied->contains($node_id)
364#    return $ways_needed->contains($node_id);
365#};
366#
367#select_ways($way_selector, $ways_copied, $nodes_needed);
368#
369#my $relation_selector = sub {
370#    my ($relation_id, $member_type, $member_id) = @_;
371#    return 0 if ($relations_copied->contains($relation_id));
372#    return 1 if ($member_type eq "node" && $nodes_copied->contains($member_id));
373#    return 1 if ($member_type eq "way" && $ways_copied->contains($member_id));
374#    return 0;
375#};
376#
377#select_relations($relation_selector, $relations_copied, $nodes_needed, $ways_needed, $relations_needed);
378
379# you may do this in a loop to fetch things recursively, using a termination
380# condition like:
381# last if ($relations_needed->subset($relations_copied) &&
382#   $ways_needed->subset($ways_copied) &&
383#   $nodes_needed->subset($nodes_copied));
384
385end_output();
386
387sub select_nodes
388{
389    my ($selector, $copied) = @_;
390    my $lastpos;
391
392    # jump to a suitable position of input
393    if (defined($filepos_nodes) && $input_is_seekable)
394    {
395        seek($input_file_handle, $filepos_nodes, 0);
396    }
397    else
398    {
399        reopen_input();
400    }
401
402    while(<$input_file_handle>) 
403    {
404        if ((/\s*<way/) && $assume_input_sorted)
405        {
406            $filepos_ways = $lastpos;
407            last;
408        }
409        last if /^\s*<\/osm>/;
410
411        if (/^\s*<osm.*version="([^"]+)"/)
412        {
413            die "this program supports 0.5 style files only; your input file has version $1"
414                unless $1 eq "0.5";
415        }
416        # Note: we allow a minus in the ID to allow incomplete files also   
417        elsif (/^\s*<node.*id=["'](-?\d+)['"].*lat=["']([Ee0-9.-]+)["'] lon=["']([Ee0-9.-]+)["']/)
418        {
419            $copy = &$selector($1, $2, $3);
420            if ($copy)
421            {
422                $copied->Bit_On($1);
423                print $output_file_handle_nodes $_;
424            }
425        }
426        elsif ($copy)
427        {
428            print $output_file_handle_nodes $_;
429        }
430        $lastpos = tell($input_file_handle) if $input_is_seekable;
431    }
432}
433
434sub select_ways
435{
436    my ($selector, $copied, $noderefs) = @_;
437    my $nodes_in_current_way = [];
438    my $copy = 0;
439    my $waybuf;
440    my $wid;
441    my $lastpos;
442
443    # jump to a suitable position of input
444    if (defined($filepos_ways) && $input_is_seekable)
445    {
446        seek($input_file_handle, $filepos_ways, 0);
447    }
448
449    while(<$input_file_handle>) 
450    {
451        if ((/\s*<relation/) && $assume_input_sorted)
452        {
453            $filepos_relations = $lastpos;
454            last;
455        }
456        last if /^\s*<\/osm>/;
457
458        if (/^\s*<way\s+id=["']([0-9-]+)['"]/)
459        {
460            $copy = 0;
461            $waybuf = $_;
462            $nodes_in_current_way = [];
463            $wid = $1;
464        }
465        elsif (/^\s*<nd ref=['"](.*)["']/)
466        {
467            if (!$copy)
468            {
469                # if way contains used node, and is not being copied yet, copy.
470                if (&$selector($wid, $1))
471                {
472                    print $output_file_handle_ways $waybuf;
473                    $copied->Bit_On($wid);
474                    $copy = 1;
475                    foreach my $nid(@$nodes_in_current_way)
476                    {
477                        $noderefs->Bit_On($nid);
478                    }
479                }
480            }
481            if ($copy)
482            {
483                print $output_file_handle_ways $_;
484                $noderefs->Bit_On($1);
485            }
486            else
487            {
488                $waybuf .= $_;
489                push(@$nodes_in_current_way, $1);
490            }
491        }
492        elsif ($copy)
493        {
494            print $output_file_handle_ways $_;
495        }
496        $lastpos = tell($input_file_handle) if $input_is_seekable;
497    }
498}
499
500sub select_relations
501{
502    my ($selector, $copied, $noderefs, $wayrefs, $relrefs) = @_;
503    my $members_in_current_relation = [];
504    my $copy = 0;
505    my $relbuf;
506    my $rid;
507    my $lastpos;
508    my $refs = { "node" => $noderefs, "way" => $wayrefs, "relation" => $relrefs };
509
510    # jump to a suitable position of input
511    if (defined($filepos_relations) && $input_is_seekable)
512    {
513        seek($input_file_handle, $filepos_relations, 0);
514    }
515
516    while(<$input_file_handle>) 
517    {
518        last if /^\s*<\/osm>/;
519
520        if (/^\s*<relation\s+id=["']([0-9-]+)['"]/)
521        {
522            $copy = 0;
523            $relbuf = $_;
524            $members_in_current_relation = [];
525            $rid = $1;
526        }
527        elsif (/^\s*<member /)
528        {
529            my ($ref, $type);
530            if (/ref=['"](.*?)['"]\s.*type=['"](.*?)['"]/)
531            {
532                ($ref, $type) = ($1, $2);
533            }
534            elsif (/type=['"](.*?)['"]\s.*ref=['"](.*?)['"]/)
535            {
536                ($ref, $type) = ($2, $1);
537            }
538            else
539            {
540                die "cannot parse member description: $_";
541            }
542
543            if (!$copy)
544            {
545                # if relation contains used object, and is not being copied yet, copy.
546                if (&$selector($rid, $type, $ref))
547                {
548                    $copied->Bit_On($rid);
549                    print $output_file_handle_rels $relbuf;
550                    $copy = 1;
551                    foreach my $mmb(@$members_in_current_relation)
552                    {
553                        my ($t, $r) = split(":", $mmb);
554                        $refs->{$t}->Bit_On($r) if defined($refs->{$t});
555                    }
556                }
557            }
558
559            if ($copy)
560            {
561                print $output_file_handle_rels $_;
562                $refs->{$type}->Bit_On($ref) if defined($refs->{$type});
563            }
564            else
565            {
566                $relbuf .= $_;
567                push(@$members_in_current_relation, "$type:$ref");
568            }
569        }
570        elsif ($copy)
571        {
572            print $output_file_handle_rels $_;
573        }
574        $lastpos = tell($input_file_handle) if $input_is_seekable;
575    }
576}
577
578sub reopen_input
579{
580    undef $input_file_handle;
581    $input_file_handle = new IO::File;
582    if ($infile =~ /\.bz2/) 
583    {
584        $input_file_handle->open("bzcat $infile |") || die "Could not bzcat file: $infile: $!";
585        $input_is_seekable = 0;
586    } 
587    elsif ($infile =~ /\.gz/) 
588    {
589        $input_file_handle->open("zcat $infile |") || die "Could not zcat file: $infile: $!";
590        $input_is_seekable = 0;
591    } 
592    else 
593    {
594        $input_file_handle->open($infile) || die "Could not open file: $infile: $!";
595        $input_is_seekable = 1;
596    }
597}
598
599sub end_output
600{
601    print $output_file_handle_rels "</osm>\n";
602    $output_file_handle_nodes->close();
603    $output_file_handle_ways->close();
604    $output_file_handle_rels->close();
605    if ($sort_output)
606    {
607        # fixme do this with seek from within perl
608        system("cat $outfile.w >> $outfile.n");
609        system("cat $outfile.r >> $outfile.n");
610        unlink("$outfile.w");
611        unlink("$outfile.r");
612        rename("$outfile.n","$outfile");
613    }
614    exit;
615}
Note: See TracBrowser for help on using the repository browser.