Ahoj,
obrázky parcel od Petra Vejsady mi nějak nedaly spát, tak jsem včera v
noci ubastlil skript na konverzi parcel do OSM.
Skriptík nemá jiné ambice než vyzkoušet si natažení parcel do JOSM a
vůbec se seznámit s VFR formátem RUIANu, takže nečekejte zázraky. Mj.
nefungují multipolygony, beru jen polygony typu gml:LinearRing apod.
Tagování jsem narychlo opsal z
http://wiki.openstreetmap.org/wiki/Cs:RUIAN, knihovny na zpracování
GML apod. jsem zatím nehledal.
Vyžaduje perl, proj, Geo::Proj4
(http://search.cpan.org/dist/Geo-Proj4/), a grid jezek_czech08.llb
přejmenovaný na "czech" (vyšťouráno tady v archívu,
http://freegis.fsv.cvut.cz/gwiki/S-JTSK_/_Grid).
Pár obrázků:
http://www.maatts.cz/obrazek/3/lomnice-ruian-parcely-png/
http://www.maatts.cz/obrazek/3/blansko-ruian-parcely-png/
http://www.maatts.cz/obrazek/3/brno-ruian-parcely-png/
Vyzkoušeno na uvedených třech obcích = na čtvrté klidně může
zhavarovat. ;-)
BACHA - ve velkých městech skript žere mraky paměti, optimalizacema
jsem se netrápil. :-)) Třeba Brno (kód obce 582786) má ve výsledku
1.122.529 uzlů a 251.328 parcel, a skriptík na to potřebuje aspoň 12GB
RAM a/nebo hodně swapu na SSD.
Martin
#!/usr/bin/perl -w
use strict;
use Geo::Proj4;
use XML::Simple;
use Data::Dumper;
#
# !!! NEPOUZIVAT PRO OSTRE MAPOVANI DO OSM, POKUSNA ALFA VERZE !!!
#
#
# Pouziti:
#
# (1) Stahnout Kompletni datovou sadu vybrane oblasti z http://vdp.cuzk.cz/vdp/ruian/vymennyformat/vyhledej, napr. podle kodu
# 581976 = Lomnice u Tisnova
# 581283 = Blansko
# 582786 = Brno (velke!)
#
# (2) gunzip 20140630_OB_581976_UKSH.xml.gz
#
# (3) par2osm.pl 20140630_OB_581976_UKSH.xml > lomnice.osm
#
#
#
my $input_file = $ARGV[0];
die "No input file specified" unless (defined $input_file);
die "Input file not found" unless (-f $input_file);
my $epsg_srs_name = 'urn:ogc:def:crs:EPSG::5514';
my $proj_from = Geo::Proj4->new('+proj=krovak +ellps=bessel +nadgrids=czech');
my $proj_to = Geo::Proj4->new('+proj=longlat +datum=WGS84');
die "Undefined source projection" unless (defined $proj_to);
die "Undefined destination projection" unless (defined $proj_from);
my $nodemap = {};
my $idgen = -1;
my $druh_pozemku_map = {
# 1 => [] ???
2 => [ 'landuse', 'farmland' ],
3 => [ 'landuse', 'hop_garden' ],
4 => [ 'landuse', 'vineyard' ],
5 => [ 'leisure', 'garden' ],
6 => [ 'landuse', 'orchard' ],
7 => [ 'landuse', 'meadow' ],
8 => [ 'landuse', 'meadow' ],
10 => [ 'landuse', 'forest', 'wood', 'yes' ],
};
my $zpusoby_vyuziti_map = {
1 => [ 'landuse', 'greenhouse_horticulture' ],
2 => [ 'landuse', 'plant_nursery' ],
3 => [ 'landuse', 'plantation', 'wood', 'yes' ],
4 => [ 'natural', 'wood', 'wood', 'yes' ],
5 => [ 'landuse', 'forest', 'wood', 'yes' ],
6 => [ 'natural', 'water', 'water', 'pond' ],
7 => [ 'waterway', 'river' ],
8 => [ 'waterway', 'canal' ],
9 => [ 'natural', 'water', 'water', 'lake' ],
10 => [ 'natural', 'water', 'water', 'reservoir' ],
11 => [ 'natural', 'wetland' ],
#12 => [ ] ???
13 => [ 'landuse', 'brownfield' ],
14 => [ 'landuse', 'railway' ],
15 => [ 'landuse', 'highway' ],
16 => [ 'landuse', 'highway' ],
17 => [ 'landuse', 'highway' ],
18 => [ 'landuse', 'transport' ],
19 => [ 'landuse', 'grass' ],
20 => [ 'landuse', 'recreation_ground' ],
21 => [ 'landuse', 'cemetery' ],
#22 => [] ???
23 => [ 'highway', 'service', 'area', 'yes' ],
24 => [ 'landuse', 'quarry' ],
25 => [ 'landuse', 'landfill' ],
#26 => [] ???
27 => [ 'natural', 'scrub' ], # ???? ... casty vyskyt
28 => [ 'natural', 'water' ],
29 => [ 'landuse', 'industrial' ],
};
sub warning
{
my ($text) = @_;
print STDERR "$text\n";
return undef;
}
sub nodemap_get
{
my ($lon, $lat) = @_;
my $key = "$lon/$lat";
my $node = $nodemap->{$key};
return $node if (defined $node);
my $krpos = [$lon, $lat];
my $osmpos = $proj_from->transform ($proj_to, $krpos);
$node = {
id => --$idgen,
krpos => $krpos,
osmpos => $osmpos,
};
$nodemap->{$key} = $node;
return $node;
}
sub parse_gml_pos
{
my ($pos) = @_;
return undef unless ($pos =~ /^\s*(\S+)\s+(\S+)\s*$/);
my $lon = $1;
my $lat = $2;
return nodemap_get ($lon, $lat);
}
sub parse_gml_poslist
{
my ($poslist) = @_;
$poslist =~ s/^\s*//;
$poslist =~ s/\s*$//;
my @list = split (/\s+/, $poslist);
die "Odd number of poslist entries: $poslist"
unless ((int(@list) % 2) == 0);
my $waynodes = [];
for (my $i = 0; $i < int(@list); $i += 2)
{
my $lon = @list[$i];
my $lat = @list[$i + 1];
push (@$waynodes, nodemap_get ($lon, $lat));
}
my $way = {
id => --$idgen,
nodes => $waynodes
};
return $way;
}
sub parse_pai_defpoint
{
my ($xdefpoint) = @_;
my $gml_point = $xdefpoint->{'gml:Point'};
return undef unless (defined $gml_point);
my $gml_pos = $gml_point->{'gml:pos'};
return undef unless (defined $gml_pos);
my $srs_name = $gml_point->{'srsName'};
return undef unless (defined $srs_name);
die "Unsupported EPSG: $srs_name" unless ($srs_name eq $epsg_srs_name);
return parse_gml_pos ($gml_pos);
}
sub parse_gml_polygon
{
my ($xpolygon, $ruian_id) = @_;
my $srs_name = $xpolygon->{'srsName'};
return undef unless (defined $srs_name);
die "Unsupported EPSG: $srs_name" unless ($srs_name eq $epsg_srs_name);
my $xexterior = $xpolygon->{'gml:exterior'};
die "$ruian_id: Polygon has no exterior"
unless (defined $xexterior);
my $xlring = $xexterior->{'gml:LinearRing'};
return warning ("$ruian_id: Polygon has no exterior linear ring")
unless (defined $xlring);
my $xposlist = $xlring->{'gml:posList'};
die "$ruian_id: Polygon has no exterior posList"
unless (defined $xposlist);
my $exterior_way = parse_gml_poslist ($xposlist);
#### TODO: parse interior ways too !!!!
my $polygon = {
id => --$idgen,
outer => [ $exterior_way ],
inner => []
};
return $polygon;
}
sub parse_parcel
{
my ($xparcel) = @_;
die "Undefined xparcel" unless (defined $xparcel);
#print Dumper $xparcel;
my $ruian_id = $xparcel->{'pai:Id'};
die "Parcel has no ID!" unless (defined $ruian_id);
my $druh_pozemku_kod = $xparcel->{'pai:DruhPozemkuKod'};
my $zpusoby_vyuziti_kod = $xparcel->{'pai:ZpusobyVyuzitiPozemku'};
my $polygon = undef;
my $defpoint = undef;
my $xgeometry = $xparcel->{'pai:Geometrie'};
if (defined $xgeometry)
{
my $xdefpoint = $xgeometry->{'pai:DefinicniBod'};
die "$ruian_id: parcel has no pai:DefinicniBod"
unless (defined $xdefpoint);
$defpoint = parse_pai_defpoint ($xdefpoint);
die "$ruian_id: cannot parse pai:DefinicniBod"
unless (defined $defpoint);
my $xorigborder = $xgeometry->{'pai:OriginalniHranice'};
if (defined $xorigborder)
{
my $xpolygon = $xorigborder->{'gml:Polygon'};
die "$ruian_id: parcel border has no gml:Polygon"
unless (defined $xpolygon);
$polygon = parse_gml_polygon ($xpolygon, $ruian_id);
}
}
my $parcel = {
id => --$idgen,
ruian_id => $ruian_id,
def_point => $defpoint,
polygon => $polygon,
druh_pozemku_kod => $druh_pozemku_kod,
zpusoby_vyuziti_kod => $zpusoby_vyuziti_kod,
};
return $parcel;
}
sub render_osm_node
{
my ($fd, $node) = @_;
return if (defined $node->{rendered});
my $node_id = $node->{id};
my $osmpos = $node->{osmpos};
print $fd " <node id=\"$node_id\" lon=\"$osmpos->[0]\" lat=\"$osmpos->[1]\">\n";
print $fd " <tag k=\"created_by\" v=\"par2osm\"/>\n";
print $fd " </node>\n";
$node->{rendered} = 1;
}
sub render_osm_way
{
my ($fd, $way, $way_tags) = @_;
my $way_id = $way->{id};
my $waynodes = $way->{nodes};
foreach my $node (@$waynodes)
{
render_osm_node ($fd, $node);
}
print $fd "<way id=\"$way_id\">\n";
print $fd " <tag k=\"created_by\" v=\"par2osm\"/>\n";
print $fd " <tag k=\"source\" v=\"cuzk:ruian\"/>\n";
print $fd " $way_tags\n";
foreach my $node (@$waynodes)
{
my $node_id = $node->{id};
print $fd " <nd ref=\"$node_id\"/>\n";
}
print $fd "</way>\n";
}
sub tarray2tags
{
my ($tarr) = @_;
my $tags = '';
for (my $i = 0; $i < int(@$tarr); $i += 2)
{
my $key = $tarr->[$i];
my $val = $tarr->[$i+1];
$tags .= "<tag k=\"$key\" v=\"$val\"/>\n";
}
return $tags;
}
sub parcel_landuse_tags
{
my ($parcel) = @_;
my $zvp = $parcel->{zpusoby_vyuziti_kod};
if (defined $zvp)
{
my $tarr = $zpusoby_vyuziti_map->{$zvp};
return tarray2tags ($tarr) if (defined $tarr);
}
my $dp = $parcel->{druh_pozemku_kod};
if (defined $dp)
{
my $tarr = $druh_pozemku_map->{$dp};
return tarray2tags ($tarr) if (defined $tarr);
}
$zvp = '?' unless (defined $zvp);
$dp = '?' unless (defined $dp);
warning ("$parcel->{ruian_id}: unknown landuse: druh_pozemku=$dp, zpusoby_vyuziti=$zvp");
#### co se zbytkem??
return
"<tag k=\"landuse\" v=\"unknown\"/>\n" .
"<tag k=\"ruian:aux:druh-pozemku\" v=\"$dp\"/>\n".
"<tag k=\"ruian:aux:zpusob-vyuziti\" v=\"$zvp\"/>\n";
}
sub render_osm_parcel
{
my ($fd, $parcel) = @_;
my $ruian_id = $parcel->{ruian_id};
my $polygon = $parcel->{polygon};
# ignore parcels without geometry...
return warning ("$ruian_id: missing geometry")
unless (defined $polygon);
my $parcel_tags = "<tag k=\"ref:ruian:parcel\" v=\"$ruian_id\"/>\n";
$parcel_tags .= parcel_landuse_tags ($parcel);
my $outers = $polygon->{outer};
my $inners = $polygon->{inner};
die "Multipolygons not supported yet!!!"
unless (int(@$outers) == 1 && int(@$inners) == 0);
my $way = $outers->[0];
render_osm_way ($fd, $way, $parcel_tags);
}
sub render_osm_xml
{
my ($fd, $parcels) = @_;
print $fd "<?xml version='1.0' encoding='UTF-8'?>\n";
print $fd "<osm version='0.6' generator='pyshp'>\n";
my $count = 0;
foreach my $parcel (@$parcels)
{
render_osm_parcel ($fd, $parcel);
++$count;
#last if ($count >= 10);
}
print $fd "<!-- parcels rendered: $count -->\n";
print $fd "</osm>\n";
}
sub main
{
my $xml = new XML::Simple;
# read XML file
my $xdata = $xml->XMLin($input_file);
my $xparcels = $xdata->{'vf:Data'}->{'vf:Parcely'}->{'vf:Parcela'};
#print Dumper $parcels;
my $parcels = [];
foreach my $xparcel (@$xparcels)
{
my $parcel = parse_parcel ($xparcel);
push (@$parcels, $parcel);
}
#### specify output file on command line?
my $fd;
open ($fd, ">&STDOUT") || die "Cannot dup stdout";
render_osm_xml ($fd, $parcels);
close ($fd);
}
main ();
1;
__END__
------- TODO ------------
* Dodelat zpracovani 'interior' polygonu a vytvaret multipolygony
* Blansko/Brno - Nezpracuji se polygony tvaru:
<gml:Polygon gml:id="HPA.19737019010" srsName="urn:ogc:def:crs:EPSG::5514" srsDimension="2">
<gml:exterior>
<gml:Ring>
<gml:curveMember>
<gml:LineString gml:id="HPA.19737019010.1">
<gml:posList>-593365.17 -1141552.76 -593367.31 -1141550.90 -593367.64 -1141547.07 -593362.82 -1141530.46 -593361.26 -1
</gml:LineString>
</gml:curveMember>
<gml:curveMember>
<gml:Curve gml:id="HPA.19737019010.2">
<gml:segments>
<gml:ArcString>
<gml:posList>-593307.73 -1141550.07 -593303.32 -1141556.05 -593301.76 -1141561.38</gml:posList>
</gml:ArcString>
</gml:segments>
</gml:Curve>
</gml:curveMember>
<gml:curveMember>
<gml:LineString gml:id="HPA.19737019010.3">
<gml:posList>-593301.76 -1141561.38 -593302.10 -1141572.54 -593304.49 -1141586.26 -593319.02 -1141583.37 -593319.60 -1
</gml:LineString>
</gml:curveMember>
</gml:Ring>
</gml:exterior>
</gml:Polygon>
_______________________________________________
Talk-cz mailing list
Talk-cz@openstreetmap.org
https://lists.openstreetmap.org/listinfo/talk-cz
