The area of a linear ring is not computed correctly

[Ari Jolma]

OGR_G_GetArea uses getGeometryType to determine the type of the geometry. For linear rings this returns linestring and thus 0 is returned.

[warmerdam]

Mateusz,

Please fix this in trunk and 1.4 branch, and add a test for this condition in gdalautotest/ogr/ogr_geom.py (trunk only is fine for the test).

[Mateusz Łoskot]

Fixed related to last reopen action have been applied (r12026, r12027 and r12028).

[Mateusz Łoskot]

I'm sorry but I've made a mistake. The comment No 3 is for Ticket #1539. I'm reopening the bug.

[Mateusz Łoskot]

Ari,

I'm testing this issue and I found not wkbLineString but some other value is returned. Could you check the exact value you're getting from getGeometryType() called on your geometry?

I created LinearRing as follows:

geom = ogr.Geometry( type = ogr.wkbLinearRing )
geom.AddPoint( 0, 0)
geom.AddPoint( 10, 0)
geom.AddPoint( 10, 10)
geom.AddPoint( 0, 10)
geom.AddPoint( 0, 0)

# Reports: -2147483646
print geom.GetGeometryType()

# Reports: 3
print geom.GetCoordinateDimension()

This type is something like a garbage. It's not even wkbLineString25D what could be expected from the coordinate dimension.

So, please provide me with some details about how do you create geometry and what exact value of geometry type you're getting. Thanks!

[Ari Jolma]

With the trunk I get:

use Geo::GDAL;

$geom = Geo::OGR::Geometry->new( $Geo::OGR::wkbLinearRing );
$geom->AddPoint( 0, 0);
$geom->AddPoint( 10, 0);
$geom->AddPoint( 10, 10);
$geom->AddPoint( 0, 10);
$geom->AddPoint( 0, 0);

print $geom->GetGeometryType(),"\n";
# Reports: 2
print $geom->GetGeometryName(),"\n";
# Reports: LINEARRING

print $geom->GetCoordinateDimension(),"\n";
# Reports: 2
print $geom->GetArea(),"\n";
# Reports: 0

2 is in fact wkbLineString since there is no GetGeometryType method in LinearRing (it inherits the one from LineString).

In the enum there is wkbLinearRing but: "non-standard, just for createGeometry()", I don't know what would be the general implications of defining GetGeometryType for LinearRing, but it would fix the GetArea problem.

[Mateusz Łoskot]

Test of geometry type and cooridnates dimensionality values for OGRLinearRing object

[Mateusz Łoskot]

Python script testing values of geometry type and cooridnates dimensionality for OGRLinearRing object

[Mateusz Łoskot]

Ari,

Your test in Perl gives correct values, except the Zero returned from GetArea().

I made two similar tests in C++ (test.cpp) and Python (ogr_linearring_type_test.py) - both tests are attached to the ticket - and there is something wrong with values I get from Python script.

• test.cpp output

mloskot@dog:~/dev/gdal/bugs/1755$ ./test 
2
2
2
100

Above, all values are correct. Please note that I fixed the OGR_G_GetArea() function using following hack (fix not applied to SVN):

...
case wkbLinearRing:
case wkbLineString:
   return ((OGRLinearRing *) hGeom)->get_Area();

• ogr_linearring_type_test.py output

mloskot@dog:~/dev/gdal/bugs/1755$ ./ogr_linearring_type_test.py 
-2147483646
3
Area: 100.000000

Above, only area value is correct, but dimension = 3 and GetGeometryType() == garbage.

I've not found where does the strange type value come from. The Python binding (pymod/gdal_wrap.cpp) looks correct to me, function _wrap_OGR_G_GetGeometryType(), it just forwards OGRGeometryH and translates enum result to integer:

   _result = (OGRwkbGeometryType )OGR_G_GetGeometryType(_arg0);

• ogr_linearring_type_test.pl - I also tested your Perl script from the comment above using GDAL with my dirty fix and it gives correctl results:

mloskot@dog:~/dev/gdal/bugs/1755$ perl ogr_linearring_type_test.pl
2
LINEARRING
2
100

The area is reported correctly, as well as coords dimension (2) and geometry type values (2).

Hmm, looks like a bug somewhere in Python bindings (gdal/pymod).

[Mateusz Łoskot]

There is something interesting and likely buggy in the gdal/pymod bindings. Check this simple Python script:

geom = ogr.Geometry( type = ogr.wkbLinearRing )
print geom.GetGeometryType()
geom.AddPoint( 0, 0)
print geom.GetGeometryType()

and here is output I'm getting:

mloskot@dog:~/dev/gdal/bugs/1755$ ./ogr_linearring_type_test.py 
2
-2147483646

The first value is correct (2 == wkbLineString), but after a point is added, value of type of geometry gets broken.

The problem is with behavior of AddPoint() method of ogr.Geometry class. This method silently sets z=0 if no z value is given. Next, it runs to call Make3D() on behalf of OGRLinearRing. After this motion, geometry type of linearring is no longer reported as wkbLineString but as wkbLineString25D, in ogrlinestring.cpp:69​. After this silent change, linearring is not properly recognized by OGR_G_GetArea() function because value of geometry type is -2147483646.

If I explicitly use AddPoint_2D(), everything looks fine:

geom = ogr.Geometry( type = ogr.wkbLinearRing )
print geom.GetGeometryType() # reports 2
geom.AddPoint_2D( 0, 0)
print geom.GetGeometryType() # reports 2 too

In Perl, the AddPoint() function behaves differently and automatically recognizes if 2 or 3D point is being added (OGR.pm:483):

@_ == 4 ? AddPoint_3D(@_) : AddPoint_2D(@_);

My question is Is it safe (by design) to add wkbLineString25D type to the switch-case in the OGR_G_GetArea() function?

[Mateusz Łoskot]

I applied simple fix to OGR_G_GetArea function (r12039) and added new test case (r12040).

[Mateusz Łoskot]

The fix has been ported to the stable branch 1.4 (r12041).