[gdal-dev] FW: Strange results of simple polygonising

Ari Jolma ari.jolma at gmail.com
Sat Jan 21 00:51:32 PST 2017


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-bounces at lists.osgeo.org] *On Behalf 
> Of *Casper Børgesen (CABO)
> *Sent:* 19. januar 2017 13:26
> *To:* gdal-dev at 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 at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170121/2bbc810e/attachment-0001.html>


More information about the gdal-dev mailing list