前言
shp转osm格式的方法有很多种，参考OpenStreetMap相关wiki：https://wiki.openstreetmap.org/wiki/Software_comparison/Import_a_shapefile
需要注意的是，国外的这些转换工具基本都要求shp文件的坐标系是WGS84的，关于如何转换坐标系我上一篇文章已经讲到过。下面介绍上图中的前两种方法完成shp文件到osm格式文件的转换。
方法一：使用Merkaartor软件导入shp后再导出osm
Merkaartor是一款地图编辑软件，支持Linux、Windows、macOS多端，可导入一种格式的地图文件然后导出多种格式，具体使用步骤如下：
1. 下载安装Merkaartor客户端
下载地址：http://merkaartor.be/p/download/
选择对应版本的安装包下载（我用的Windows 64位），然后直接安装就行了
2. 导入shp文件然后导出osm格式
安装完成后打开Merkaartor客户端，选择File->Import using GDAL，在弹出的窗口中点Yes，然后选择自己要转换的shp文件，等待文件导入
导入完成后在左侧点击选中自己的shp文件，然后点击选择File->Export->OSM，进行导出即可。
但是这种方法的缺点就是Merkaartor会自动给属性字段名称加上两个下划线，这点有时候就很蛋疼。转换出的格式如下
方法二：使用perl脚本shp2osm对文件进行转换
这里我们使用shp2osm对地图进行转换。shp2osm是一个perl脚本，文档地址：https://wiki.openstreetmap.org/wiki/Shp2osm
这种方法可以非常快速便捷的实现shp到osm格式的转换，但是这是个2008年写的脚本，其转换出的osm文件的格式和现在版本的osm格式可能还是稍微有些区别，虽然也能在Arcgis中打开，但是有些操作就会产生不兼容的情况（比如尝试用Osmosis将它导入到Postgresql数据库中），有不适配的地方就只能自己修改脚本
具体使用步骤如下：
1. 下载shp2osm脚本并保存
确保shp文件的坐标系是WGS84后，在shp文件所在的目录下建立一个shp2osm.pl文件，然后将以下源码全部复制进去然后保存。
shp2osm.pl源码地址：https://svn.openstreetmap.org/applications/utils/import/shp2osm/shp2osm.pl
为了防丢失我直接复制一份在这里吧
# Copyright (c) 2006 Gabriel Ebner <ge@gabrielebner.at>
# updated in 2008 by Tobias Wendorff <tobias.wendorff@uni-dortmund.de>
# HTML-Entities based on ideas of Hermann Schwärzler
# Gauß-Krüger implementation based on gauss.pl by Andreas Achtzehn
# version 1.3 (17. September 2008)
use Geo::ShapeFile;
use HTML::Entities qw(encode_entities_numeric);
use Math::Trig;
if(@ARGV == 0) {
print "usage:\n";
print "with transformation from GK to WGS84: 'shp2osm.pl shapefile GK'\n";
print "without transformation: 'shp2osm.pl shapefile'";
exit;
}
print <<'END';
<?xml version='1.0'?>
<osm version='0.5' generator='shp2osm.pl'>
END
#BEGIN { our %default_tags = ( natural => 'water', source => 'SWDB' ); }
BEGIN { our %default_tags = ( ); }
my $file = @ARGV[0];
$file =~ s/\.shp$//;
my $shpf = Geo::ShapeFile->new($file);
proc_shpf($shpf);
{
BEGIN { our $i = -1; }
sub tags_out {
my ($tags) = @_;
my %tags = $tags ? %$tags : ();
#$tags{'created_by'} ||= 'shp2osm.pl';
delete $tags{'_deleted'} unless $tags{'_deleted'};
while ( my ( $k, $v ) = each %tags ) {
my $key = encode_entities_numeric($k);
my $val = encode_entities_numeric($v);
print ' <tag k="'. $key .'" v="'. $val ."\"/>\n" if $val;
}
}
sub node_out {
my ( $lon, $lat, $tags ) = @_;
my $id = $i--;
if(@ARGV[1] eq 'GK') {
my ($wgs84lon, $wgs84lat) = gk2geo($lon, $lat);
print " <node id='$id' visible='true' lat='$wgs84lat' lon='$wgs84lon' />\n";
} else {
print " <node id='$id' visible='true' lat='$lat' lon='$lon' />\n";
}
$id;
}
sub seg_out {
my $id = $i+1;
$id;
}
sub way_out {
my ( $segs, $tags ) = @_;
my $id = $i--;
print " <way id='$id' visible='true'>\n";
print " <nd ref='$_' />\n" for @$segs;
tags_out $tags;
print " </way>\n";
$id;
}
}
print <<'END';
</osm>
END
sub polyline_out {
my ( $pts, $tags, $connect_last_seg ) = @_;
my ( $first_node, $last_node, @segs );
for my $pt (@$pts) {
my $node = node_out @$pt;
push @segs, seg_out $last_node, $node;
$last_node = $node;
$first_node ||= $last_node;
}
push @segs, seg_out $last_node, $first_node
if $first_node && $connect_last_seg;
way_out \@segs, $tags;
}
sub proc_obj {
my ( $shp, $dbf, $type ) = @_;
my $tags = { %default_tags, %$dbf };
my $is_polygon = $type % 10 == 5;
for ( 1 .. $shp->num_parts ) {
polyline_out [ map( [ $_->X(), $_->Y() ], $shp->get_part($_) ) ], $tags,
$is_polygon;
}
}
sub proc_shpf {
my ($shpf) = @_;
my $type = $shpf->shape_type;
for ( 1 .. $shpf->shapes() ) {
my $shp = $shpf->get_shp_record($_);
my %dbf = $shpf->get_dbf_record($_);
proc_obj $shp, \%dbf, $type;
}
}
sub gk2geo {
my $sr = $_[0];
my $sx = $_[1];
my $bm = int($sr/1000000);
my $y = $sr-($bm*1000000+500000);
my $si = $sx/111120.6196;
my $px = $si+0.143885358*sin(2*$si*0.017453292)+0.00021079*sin(4*$si*0.017453292)+0.000000423*sin(6*$si*0.017453292);
my $t = (sin($px*0.017453292))/(cos($px*0.017453292));
my $v = sqrt(1+0.006719219*cos($px*0.017453292)*cos($px*0.017453292));
my $ys = ($y*$v)/6398786.85;
my $lat = $px-$ys*$ys*57.29577*$t*$v*$v*(0.5-$ys*$ys*(4.97-3*$t*$t)/24);
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)));
my $lon = $bm*3+$dl;
my $potsd_a = 6377397.155;
my $wgs84_a = 6378137.0;
my $potsd_f = 1/299.152812838;
my $wgs84_f = 1/298.257223563;
my $potsd_es = 2*$potsd_f - $potsd_f*$potsd_f;
my $potsd_dx = 606.0;
my $potsd_dy = 23.0;
my $potsd_dz = 413.0;
my $pi = 3.141592654;
my $latr = $lat/180*$pi;
my $lonr = $lon/180*$pi;
my $sa = sin($latr);
my $ca = cos($latr);
my $so = sin($lonr);
my $co = cos($lonr);
my $bda = 1-$potsd_f;
my $delta_a = $wgs84_a - $potsd_a;
my $delta_f = $wgs84_f - $potsd_f;
my $rn = $potsd_a / sqrt(1-$potsd_es*sin($latr)*sin($latr));
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)));
my $ta = (-$potsd_dx*$sa*$co - $potsd_dy*$sa*$so)+$potsd_dz*$ca;
my $tb = $delta_a*(($rn*$potsd_es*$sa*$ca)/$potsd_a);
my $tc = $delta_f*($rm/$bda+$rn*$bda)*$sa*$ca;
my $dlat = ($ta+$tb+$tc)/$rm;
my $dlon = (-$potsd_dx*$so + $potsd_dy*$co)/($rn*$ca);
my $wgs84lat = ($latr + $dlat)*180/$pi;
my $wgs84lon = ($lonr + $dlon)*180/$pi;
return( $wgs84lon, $wgs84lat);
}
2. 安装perl和相关扩展包
shp2osm.pl是perl程序所以肯定要安装perl才能运行啦，perl安装包下载地址：
https://strawberryperl.com/
安装完成再在cmd中用cpan命令安装一个操作shp文件的Perl扩展包
cpan Geo::ShapeFile
3. 运行shp2osm.pl生成osm格式的文件
打开cmd（不能用PowerShell，只有cmd中才能执行下面的DOS命令），定位到shp文件所在的目录（之前已经在这里面生成了shp2osm.pl文件），执行如下命令（*.shp即表示匹配该目录下所有shp文件）：
FOR /R .\ %G IN (*.shp) DO shp2osm.pl "%~dpnG" > "%~dpnG.osm"
等待转换完成
4. 字段的转换问题
在Arcgis中打开转换出的osm文件查看道路的属性表可能会发现，之前shp文件里的字段这里都没有，只有道路的形状信息，这是因为这里只显示合法的字段（也就是字段名称与osm格式一致的字段）,
下图是我直接从转换出的osm文件中截出的一部分数据，如图所示，osm就是用xml的格式存储的道路数据，其中一个way节点就存储着一条道路的信息，其中一条tag就是道路的一个字段属性值，shp文件中的字段都转换成了这样的tag，但是名称与osm格式中不一致的都是非法tag，所以如果有需要的话可以先将shp文件中的字段名称与osm格式也保持一致，这样最后在Arcgis的属性表中就能看到了。关于如何修改shp文件中的字段名称可参考：https://jingyan.baidu.com/article/4f34706e9ca2d7a387b56dd5.html
2810
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
