Casper,
The polygonizer can be used in 4 or 8 connected modes. I'm not sure
which one you used. My results are the following. (I used GeoTransform
(0,1,0,3,0,-1) for the rasters.
4-connected:
1:
POLYGON ((0 3,0 0,2 0,2 1,1 1,1 2,2 2,2 1,3 1,3 3,0 3))
POLYGON ((1 2,1 1,2 1,2 2,1 2))
POLYGON ((2 1,2 0,3 0,3 1,2 1))
2:
POLYGON ((0 3,0 0,3 0,3 2,2 2,2 1,1 1,1 2,2 2,2 3,0 3))
POLYGON ((2 3,2 2,3 2,3 3,2 3))
POLYGON ((1 2,1 1,2 1,2 2,1 2))
3:
POLYGON ((0 3,0 2,1 2,1 3,0 3))
POLYGON ((1 3,1 2,0 2,0 0,2 0,3 0,3 3,1 3),(1 2,1 1,2 1,2 2,1 2))
POLYGON ((1 2,1 1,2 1,2 2,1 2))
4:
POLYGON ((1 2,1 1,2 1,2 2,1 2))
POLYGON ((0 1,0 0,1 0,1 1,0 1))
POLYGON ((0 3,0 1,1 1,1 0,3 0,3 3,0 3),(1 2,2 2,2 1,1 1,1 2))
3+4:
POLYGON ((0 4,0 3,1 3,1 4,0 4))
POLYGON ((0 3,0 1,1 1,1 3,0 3))
POLYGON ((1 3,1 1,2 1,2 3,1 3))
POLYGON ((0 1,0 0,1 0,1 1,0 1))
POLYGON ((1 4,1 3,2 3,2 1,1 1,1 0,2 0,3 0,3 4,1 4))
8-connected:
1:
POLYGON ((0 3,0 0,2 0,2 1,1 1,1 2,2 2,2 1,3 1,3 3,0 3))
POLYGON ((1 2,1 1,2 1,2 0,3 0,3 1,2 1,2 2,1 2))
2:
POLYGON ((0 3,0 0,3 0,3 2,2 2,2 1,1 1,1 2,2 2,2 3,0 3))
POLYGON ((2 3,2 2,1 2,1 1,2 1,2 2,3 2,3 3,2 3))
3:
POLYGON ((0 3,0 2,1 2,1 1,2 1,2 2,1 2,1 3,0 3))
POLYGON ((1 3,1 2,0 2,0 0,2 0,3 0,3 3,1 3),(1 2,1 1,2 1,2 2,1 2))
4:
POLYGON ((0 3,0 1,1 1,1 0,3 0,3 3,0 3),(1 2,2 2,2 1,1 1,1 2))
POLYGON ((1 2,1 1,0 1,0 0,1 0,1 1,2 1,2 2,1 2))
3+4:
POLYGON ((0 4,0 3,1 3,1 1,0 1,0 0,1 0,1 1,2 1,2 3,1 3,1 4,0 4))
POLYGON ((1 4,1 3,0 3,0 1,1 1,1 0,2 0,3 0,3 4,1 4),(1 3,1 1,2 1,2 3,1 3))
To me the results seem to be correct by manual checking with QGIS.
Can you show the code that you used?
Best,
Ari
the code that I used to generate these results:
use Modern::Perl;
use Geo::GDAL;
sub polygonize {
my ($b, $name) = @_;
my $l2 = Geo::OGR::Driver('ESRI
Shapefile')->Create('.')->CreateLayer($name);
my $options = {'8CONNECTED' => 1};
my $l = $b->Polygonize(Options => $options);
$l->ResetReading();
while (my $f = $l->GetNextFeature()) {
$l2->CreateFeature($f);
my $g = $f->GetGeometry();
say $g->As(Format => 'WKT');
}
}
system "rm test*";
my $b = Geo::GDAL::Driver('MEM')->Create(Width => 3, Height => 3)->Band;
$b->Dataset->GeoTransform(0,1,0,3,0,-1);
$b->Dataset->SpatialReference(Geo::OSR::SpatialReference->new(EPSG=>3067));
my $d = $b->ReadTile;
# 1:
$d->[0] = [1,1,1];
$d->[1] = [1,0,1];
$d->[2] = [1,1,0];
$b->WriteTile($d);
$b->Dataset->Translate('test1.tiff');
say "1:";
polygonize($b, 'test1');
$d->[0] = [1,1,0];
$d->[1] = [1,0,1];
$d->[2] = [1,1,1];
$b->WriteTile($d);
$b->Dataset->Translate('test2.tiff');
say "2:";
polygonize($b, 'test2');
$d->[0] = [0,1,1];
$d->[1] = [1,0,1];
$d->[2] = [1,1,1];
$b->WriteTile($d);
$b->Dataset->Translate('test3.tiff');
say "3:";
polygonize($b, 'test3');
$d->[0] = [1,1,1];
$d->[1] = [1,0,1];
$d->[2] = [0,1,1];
$b->WriteTile($d);
$b->Dataset->Translate('test4.tiff');
say "4:";
polygonize($b, 'test4');
$b = Geo::GDAL::Driver('MEM')->Create(Width => 3, Height => 4)->Band;
$b->Dataset->GeoTransform(0,1,0,4,0,-1);
$b->Dataset->SpatialReference(Geo::OSR::SpatialReference->new(EPSG=>3067));
$d = $b->ReadTile;
$d->[0] = [0,1,1];
$d->[1] = [1,0,1];
$d->[2] = [1,0,1];
$d->[3] = [0,1,1];
$b->WriteTile($d);
$b->Dataset->Translate('test3+4.tiff');
say "3+4:";
polygonize($b, 'test3+4');
20.01.2017, 10:31, Casper Børgesen (CABO) kirjoitti:
I have looked at the special case of a raster combining situation 3
and 4 below:
3 + 4:
# #
# #
# #
# #
The polygoniser handles this fine by return a result with two separate
polygons.
For easier reference, I have included the WKT of the 4 + 1 situations
and their results:
1:
Input: POLYGON ((1 1,3 1,3 2,2 2,2 3,3 3,3 2,4 2,4 4,1 4,1 1))
Result: POLYGON ((1 4,1 1,3 1,3 2,2 2,2 3,3 3,3 2,4 2,4 4,1 4))
2:
Input: POLYGON ((1 1,4 1,4 3,3 3,3 2,2 2,2 3,3 3,3 4,1 4,1 1))
Result: POLYGON ((1 4,1 1,4 1,4 3,3 3,3 2,2 2,2 3,3 3,3 4,1 4))
3:
Input: POLYGON ((1 1,4 1,4 4,2 4,2 3,3 3,3 2,2 2,2 3,1 3,1 1))
Result: POLYGON ((2 4,2 3,1 3,1 1,3 1,4 1,4 4,2 4),(2 3,2 2,3 2,3 3,2 3))
4:
Input: POLYGON ((1 4,1 2,2 2,2 3,3 3,3 2,2 2,2 1,4 1,4 4,1 4))
Result: POLYGON ((1 4,1 2,2 2,2 1,4 1,4 4,1 4),(2 3,3 3,3 2,2 2,2 3))
3 + 4:
Input: MULTIPOLYGON (((1 2,2 2,2 4,1 4,1 2)),((2 1,4 1,4 5,2 5,2 4,3
4,3 2,2 2,2 1)))
Result: POLYGON ((1 4,1 2,2 2,2 4,1 4)) + POLYGON ((2 5,2 4,3 4,3 2,2
2,2 1,3 1,4 1,4 5,2 5))
I have performed the tests on both GDAL 1.11 and GDAL 2.1 and both
return the same results.
Can anyone confirm my observations?
/Background: The 4 situations are very simplified versions of much
more complicated polygons, where the problem occurs. If I save the
results from situation 3 or 4 to a shape file with ogr, the hole in
the results are converted to an exterior ring in a separate polygon
thus creating overlapping results./
Regards, Casper
*From:*gdal-dev [mailto:gdal-dev-boun...@lists.osgeo.org] *On Behalf
Of *Casper Børgesen (CABO)
*Sent:* 19. januar 2017 13:26
*To:* gdal-dev@lists.osgeo.org
*Subject:* [gdal-dev] Strange results of simple polygonising
Hi,
I have four simple rasters that I would like to polygonise:
1:
# # #
# #
# #
2:
# #
# #
# # #
3:
# #
# #
# # #
4:
# # #
# #
# #
It’s the same shape, just rotated 90, 180, 270 degrees. The resulting
polygons for 1 and 2 are simple polygons without holes, where the
results for 3 and 4 are polygons with a hole intersecting the exterior
ring. Thus 3 and 4 results in invalid geometries with self-intersection.
Is this the intended behavior of the polygoniser?
Regards, Casper
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
