source: subversion/applications/utils/import/shp2osm/shp2osm.pl @ 20034

Last change on this file since 20034 was 10856, checked in by tobwen, 12 years ago

Shapefile to OSM converter

File size: 4.8 KB
Line 
1# Copyright (c) 2006 Gabriel Ebner <ge@gabrielebner.at>
2# updated in 2008 by Tobias Wendorff <tobias.wendorff@uni-dortmund.de>
3# HTML-Entities based on ideas of Hermann Schwärzler
4# Gauß-Krüger implementation based on gauss.pl by Andreas Achtzehn
5# version 1.3 (17. September 2008)
6
7use Geo::ShapeFile;
8use HTML::Entities qw(encode_entities_numeric);
9use Math::Trig;
10
11if(@ARGV == 0) {
12   print "usage:\n";
13   print "with transformation from GK to WGS84: 'shp2osm.pl shapefile GK'\n";
14   print "without transformation: 'shp2osm.pl shapefile'";
15   exit;
16}
17
18print <<'END';
19<?xml version='1.0'?>
20<osm version='0.5' generator='shp2osm.pl'>
21END
22
23#BEGIN { our %default_tags = ( natural => 'water', source => 'SWDB' ); }
24BEGIN { our %default_tags = (  ); }
25
26my $file = @ARGV[0];
27$file =~ s/\.shp$//;
28my $shpf = Geo::ShapeFile->new($file);
29proc_shpf($shpf);
30
31{
32     BEGIN { our $i = -1; }
33
34     sub tags_out {
35         my ($tags) = @_;
36         my %tags = $tags ? %$tags : ();
37         #$tags{'created_by'} ||= 'shp2osm.pl';
38         delete $tags{'_deleted'} unless $tags{'_deleted'};
39         while ( my ( $k, $v ) = each %tags ) {
40             my $key = encode_entities_numeric($k);
41             my $val = encode_entities_numeric($v);
42             print '    <tag k="'. $key .'" v="'. $val ."\"/>\n" if $val;
43         }
44     }
45
46     sub node_out {
47         my ( $lon, $lat, $tags ) = @_;
48         my $id = $i--;
49         if(@ARGV[1] eq 'GK') {
50             my ($wgs84lon, $wgs84lat) = gk2geo($lon, $lat);
51             print "  <node id='$id' visible='true' lat='$wgs84lat' lon='$wgs84lon' />\n";
52         } else {
53            print "  <node id='$id' visible='true' lat='$lat' lon='$lon' />\n";         
54         }
55         $id;
56     }
57
58     sub seg_out {
59         my $id = $i+1;
60         $id;
61     }
62
63     sub way_out {
64         my ( $segs, $tags ) = @_;
65         my $id = $i--;
66         print "  <way id='$id' visible='true'>\n";
67         print "    <nd ref='$_' />\n" for @$segs;
68         tags_out $tags;
69         print "  </way>\n";
70         $id;
71     }
72}
73
74
75print <<'END';
76</osm>
77END
78
79sub polyline_out {
80    my ( $pts, $tags, $connect_last_seg ) = @_;
81
82    my ( $first_node, $last_node, @segs );
83    for my $pt (@$pts) {
84        my $node = node_out @$pt;
85        push @segs, seg_out $last_node, $node;
86        $last_node = $node;
87        $first_node ||= $last_node;
88    }
89    push @segs, seg_out $last_node, $first_node
90      if $first_node && $connect_last_seg;
91    way_out \@segs, $tags;
92}
93
94sub proc_obj {
95    my ( $shp, $dbf, $type ) = @_;
96    my $tags = { %default_tags, %$dbf };
97    my $is_polygon = $type % 10 == 5;
98    for ( 1 .. $shp->num_parts ) {
99        polyline_out [ map( [ $_->X(), $_->Y() ], $shp->get_part($_) ) ], $tags,
100          $is_polygon;
101    }
102 }
103
104sub proc_shpf {
105    my ($shpf) = @_;
106    my $type = $shpf->shape_type;
107    for ( 1 .. $shpf->shapes() ) {
108        my $shp = $shpf->get_shp_record($_);
109        my %dbf = $shpf->get_dbf_record($_);
110        proc_obj $shp, \%dbf, $type;
111    }
112}
113
114sub gk2geo {
115  my $sr  = $_[0];
116  my $sx  = $_[1];
117
118  my $bm  = int($sr/1000000);
119  my $y   = $sr-($bm*1000000+500000);
120  my $si  = $sx/111120.6196;
121  my $px  = $si+0.143885358*sin(2*$si*0.017453292)+0.00021079*sin(4*$si*0.017453292)+0.000000423*sin(6*$si*0.017453292);
122  my $t   = (sin($px*0.017453292))/(cos($px*0.017453292));
123  my $v   = sqrt(1+0.006719219*cos($px*0.017453292)*cos($px*0.017453292));
124  my $ys  = ($y*$v)/6398786.85;
125  my $lat = $px-$ys*$ys*57.29577*$t*$v*$v*(0.5-$ys*$ys*(4.97-3*$t*$t)/24);
126  my $dl  = $ys*57.29577/cos($px*0.017453292) * (1-$ys*$ys/6*($v*$v+2*$t*$t-$ys*$ys*(0.6+1.1*$t*$t)*(0.6+1.1*$t*$t)));
127  my $lon = $bm*3+$dl;
128
129  my $potsd_a  = 6377397.155;
130  my $wgs84_a  = 6378137.0;
131  my $potsd_f  = 1/299.152812838;
132  my $wgs84_f  = 1/298.257223563;
133
134  my $potsd_es = 2*$potsd_f - $potsd_f*$potsd_f;
135
136  my $potsd_dx = 606.0;
137  my $potsd_dy = 23.0;
138  my $potsd_dz = 413.0;
139  my $pi = 3.141592654;
140  my $latr = $lat/180*$pi;
141  my $lonr = $lon/180*$pi;
142
143  my $sa = sin($latr);
144  my $ca = cos($latr);
145  my $so = sin($lonr);
146  my $co = cos($lonr);
147
148  my $bda  = 1-$potsd_f;
149
150  my $delta_a = $wgs84_a - $potsd_a;
151  my $delta_f = $wgs84_f - $potsd_f;
152
153  my $rn = $potsd_a / sqrt(1-$potsd_es*sin($latr)*sin($latr));
154  my $rm = $potsd_a * ((1-$potsd_es)/sqrt(1-$potsd_es*sin($latr)*sin($latr)*1-$potsd_es*sin($latr)*sin($latr)*1-$potsd_es*sin($latr)*sin($latr)));
155
156  my $ta = (-$potsd_dx*$sa*$co - $potsd_dy*$sa*$so)+$potsd_dz*$ca;
157  my $tb = $delta_a*(($rn*$potsd_es*$sa*$ca)/$potsd_a);
158  my $tc = $delta_f*($rm/$bda+$rn*$bda)*$sa*$ca;
159  my $dlat = ($ta+$tb+$tc)/$rm;
160
161  my $dlon = (-$potsd_dx*$so + $potsd_dy*$co)/($rn*$ca);
162
163  my $wgs84lat = ($latr + $dlat)*180/$pi;
164  my $wgs84lon = ($lonr + $dlon)*180/$pi;
165
166  return( $wgs84lon, $wgs84lat);
167}
Note: See TracBrowser for help on using the repository browser.