Ticket #1755 (closed defect: fixed)

Opened 6 years ago

Last modified 6 years ago

The area of a linear ring is not computed correctly

Reported by: ajolma Owned by: mloskot
Priority: high Milestone: 1.4.3
Component: OGR_SF Version: unspecified
Severity: normal Keywords: geometry area
Cc: warmerdam

Description

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

Attachments

test.cpp Download (0.6 KB) - added by mloskot 6 years ago.
Test of geometry type and cooridnates dimensionality values for OGRLinearRing object
ogr_linearring_type_test.py Download (321 bytes) - added by mloskot 6 years ago.
Python script testing values of geometry type and cooridnates dimensionality for OGRLinearRing object

Change History

Changed 6 years ago by warmerdam

  • cc warmerdam added
  • keywords geometry area added
  • owner changed from warmerdam to mloskot
  • milestone set to 1.4.3

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).

Changed 6 years ago by mloskot

  • priority changed from normal to high
  • status changed from new to assigned

Changed 6 years ago by mloskot

  • status changed from assigned to closed
  • resolution set to fixed

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

Changed 6 years ago by mloskot

  • status changed from closed to reopened
  • resolution fixed deleted

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

Changed 6 years ago by mloskot

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!

Changed 6 years ago by ajolma

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.

Changed 6 years ago by mloskot

Test of geometry type and cooridnates dimensionality values for OGRLinearRing object

Changed 6 years ago by mloskot

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

Changed 6 years ago by mloskot

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).

Changed 6 years ago by mloskot

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?

Changed 6 years ago by mloskot

  • status changed from reopened to closed
  • resolution set to fixed

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

Changed 6 years ago by mloskot

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

Note: See TracTickets for help on using tickets.