Bonjour Le site du cadastre permet d'exporter une zone au format pdf. En réalisant cet export sur une bounding box suffisamment grande, on peut obtenir la carte entière d'une ville dans un format vectoriel.
Avec un script perl un peu gruik, j'ai pu extraire d'un coup le bâti d'une planche cadastral. Les relations multipolygon sont gérées, un exemple avec le 1er arrondissement de Paris : http://dl.free.fr/foJJrx89i Ou carrément tous les arrondissements (250Mo décompressé) : http://dl.free.fr/mpAGvHFhw Je joins le script, j'utilise pdf2svg et les modules perl XML::Parser et Geo::OSR.
import-bati.sh
Description: Bourne shell script
#!/usr/bin/perl
# Usage: svg-parser.pl [EPSG] [fichier.svg] [bbox]
use XML::Parser;
use Geo::OSR;
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
my $tag_source = "cadastre-dgi-fr source : Direction Générale des Impôts - Cadastre. Mise à jour : " . ($year + 1900);
my @bbox_lbr93;
my @bbox_pts;
# Hash qui à chaque point associe une ref
my %points;
my @ways;
my @relations;
my $refnodes = 0;
my $refways = 0;
my $refrel = 0;
my $source = Geo::OSR::SpatialReference->create ();
my $target = Geo::OSR::SpatialReference->create ();
$source->ImportFromEPSG ($ARGV[0]);
$target->ImportFromEPSG ('4326');
my $transf = new Geo::OSR::CoordinateTransformation ($source, $target);
my $parser = new XML::Parser ( Handlers => {
Start => \&hdl_start,
End => \&hdl_end,
Default => \&hdl_def
});
my $surface = 0;
sub hdl_start {
my ($p, $elt, %atts) = @_;
@bbox_lbr93 = ($ARGV[2],$ARGV[3],$ARGV[4],$ARGV[5]);
$surface = 1 if ($surface == 0 && $elt eq 'g' && $atts{'id'} eq 'surface0');
if ($surface && $elt eq 'path' )
{
my @m = get_matrix ($atts{'transform'});
if ($atts{'style'} =~ m/fill:rgb\(100%,[^0]/)
{
my $s = $atts{'d'};
do {
my @points = format_point (\$s,@m);
if (!defined(@bbox_pts) && ($atts{'style'} =~ m/rgb\(100%,100%,100%\)/))
{
@bbox_pts = minmax(@points);
}
elsif (defined(@bbox_pts) && defined(@points))
{
$relations[$refrel] .= "$refways ";
new_way(@points);
}
} while ($s =~ m/M (-?\d*\.?\d*) (-?\d*\.?\d*) L/);
$refrel++;
}
}
}
sub hdl_end {}
sub hdl_def {}
sub format_point {
my ($s,@m) = @_;
my @points;
return unless $$s =~ s/^M //;
$points[0] = transform_point ($s,@m);
my $i = 1;
while ($$s =~ s/^L //)
{
$points[$i] = transform_point ($s,@m);
$i += 1;
}
$$s =~ s/^Z //;
return @points
}
sub transform_point {
my ($s,@m) = @_;
my $p;
($p->[0],$p->[1]) = $$s =~ m/(-?\d*\.?\d*) (-?\d*\.?\d*) ?/;
$$s =~ s/(-?\d*\.?\d*) (-?\d*\.?\d*) ?//;
($p->[0],$p->[1]) = (($m[0]*$p->[0] + $m[2]*$p->[1] + $m[4]),($m[1]*$p->[0] + $m[3]*$p->[1] + $m[5]));
if (defined(@bbox_pts))
{
($p->[0],$p->[1]) = (
(($p->[0] - $bbox_pts[0]) * ($bbox_lbr93[2]-$bbox_lbr93[0])/($bbox_pts[2]-$bbox_pts[0]) + $bbox_lbr93[0]),
(($p->[1] - $bbox_pts[1]) * ($bbox_lbr93[1]-$bbox_lbr93[3])/($bbox_pts[3]-$bbox_pts[1]) + $bbox_lbr93[3])
);
return ($transf->TransformPoint ($p->[0],$p->[1]));
}
else {
return $p;
}
}
sub get_matrix {
my ($s) = @_;
return split (/,/, $s) if $s =~ s/matrix\((.*)\)/$1/;
return (1,0,0,1,0,0)
}
sub minmax {
my @nodes = @_;
my $xmin;
my $ymin;
my $xmax;
my $ymax;
my $node;
foreach $node (@nodes) {
$xmin = $node->[0] if ($node->[0] < $xmin);
$ymin = $node->[1] if ($node->[1] < $ymin);
$xmax = $node->[0] if ($node->[0] > $xmax);
$ymax = $node->[1] if ($node->[1] > $ymax);
}
return ($xmin,$ymin,$xmax,$ymax)
}
sub new_point {
my ($lat,$lon) = @_;
if (!defined($points{"$lon,$lat"}))
{
$points{"$lon,$lat"} = $refnodes;
$refnodes++;
}
$ref = $points{"$lon,$lat"};
}
sub new_way {
my (@points) = @_;
foreach $node (@points) {
($lon,$lat) = ($node->[0],$node->[1]);
new_point($lat,$lon);
my $ref = $points{"$lon,$lat"};
$ways[$refways] .= "$ref ";
}
$refways++;
}
$parser->parsefile($ARGV[1]);
my $bbox_wgs84_min = ($transf->TransformPoint (@bbox_lbr93[0,1]));
my $bbox_wgs84_max = ($transf->TransformPoint (@bbox_lbr93[2,3]));
print "<?xml version='1.0' encoding='UTF-8'?>\n";
print "<osm version='0.6' generator='plop'>\n";
print "<bounds minlat=\"$bbox_wgs84_min->[1]\" minlon=\"$bbox_wgs84_min->[0]\" maxlat=\"$bbox_wgs84_max->[1]\" maxlon=\"$bbox_wgs84_max->[0]\"/>\n";
foreach $node (keys %points)
{
($lon,$lat) = split (/,/, $node);
$ref = -1 - $points{$node};
print "<node id=\"$ref\" lat=\"$lat\" lon=\"$lon\">\n";
print " <tag k=\"source\" v=\"$tag_source\"/>\n";
print "</node>\n";
}
foreach $i (0..$#relations)
{
my @relways = (split (/ /,$relations[$i]));
if (@relways >1) {
foreach $i (@relways)
{
my $refway = -1-$i;
print "<way id=\"$refway\">\n";
foreach $node (split (/ /,$ways[$i]))
{
my $ref = -1 - $node;
print " <nd ref=\"$ref\"/>\n";
}
print " <tag k=\"source\" v=\"$tag_source\"/>\n";
print " <tag k=\"building\" v=\"yes\"/>\n" if ($i == $relways[0]);
print "</way>\n";
}
my $refrel = -1-$i;
print "<relation id=\"$refrel\">\n";
print " <tag k=\"type\" v=\"multipolygon\"/>\n";
foreach $way (@relways)
{
my $ref = -1-$way;
if ($way == $relways[0]) {
print " <member type=\"way\" ref=\"$ref\" role=\"outer\"/>\n";
}
else {
print " <member type=\"way\" ref=\"$ref\" role=\"inner\"/>\n";
}
}
print "</relation>\n";
}
else
{
my $refway = -1 - $relways[0];
print "<way id=\"$refway\">\n";
foreach $node (split (/ /,$ways[$relways[0]]))
{
my $ref = -1 - $node;
print " <nd ref=\"$ref\"/>\n";
}
print " <tag k=\"source\" v=\"$tag_source\"/>\n";
print " <tag k=\"building\" v=\"yes\"/>\n";
print "</way>\n";
}
}
print "</osm>";
_______________________________________________ Talk-fr mailing list Talk-fr@openstreetmap.org http://lists.openstreetmap.org/listinfo/talk-fr
