[gdal-dev] Re: OGR/ GEOS GetNextFeature problem

[Frank Warmerdam]

Ahmet Murat Özdemiray wrote:
> Hi,
> I have a problem with OGRLayer.GetNextFeature() method. I did not
> where to post this, so here I am.

Murat,

This would be a good topic to discuss on the gdal-dev mailing
list.  I am taking the liberty of cc:ing my reply there since the
issue does not seem to be confidential in any way.

> I had a working display program using gdal, which reads ESRI shape
> files and displays contents of it on a Windows Form. I downloaded gdal
> version 1.4.3 and compiled on Win32 using nmake. It was working fine,
> until I realized that, geometric operations like Contains, Disjoints,
> Intersects etc. does not work without GEOS library. So I downloaded
> lastest versions of gdal (1.5.0) and geos(3.0.0), compiled both of
> them with release configuration for Win32 with VS2005.
> 
> When I updated my old program, I realized that things are displayed
> slower. Started debugging and saw that GetNextFeature takes up to 15
> secs for some features and a shape file containing all countries of
> the world takes 2 minutes to display. Then, I downloaded FWTools
> 2.0.4, changed all bindings to FWTools' libraries and I got similar
> results. However, when I try to open same file with OpenEV, it can
> display it in 2 secs.
> 
> Here's are a few lines of my code, what do I do wrong?
> 
> layer->SetSpatialFilterRect(-180, -90, 180, 90);
> layer->ResetFeature();

I assume the above is actually layer->ResetReading()?

> OGRFeature* ftr;
> while((ftr = layer->GetNextFeature()) != NULL)
> {
>    // do something
>    OGRFeature::DestroyFeature( ftr );
> }

You haven't done anything wrong, but you have learned that GEOS
intersection tests can be inordinately expensive.

Hmm, on further reflection reading from a shapefile using a rectangular
filter such as any filter set with SetSpatialFilterRect() should short
circuit use of GEOS as per the code in OGRLayer::FilterGeometry():

/* -------------------------------------------------------------------- */
/*      Fallback to full intersect test (using GEOS) if we still        */
/*      don't know for sure.                                            */
/* -------------------------------------------------------------------- */
     if( m_bFilterIsEnvelope )
         return TRUE;
     else
     {
         if( OGRGeometryFactory::haveGEOS() )
             return m_poFilterGeom->Intersects( poGeometry );
         else
             return TRUE;
     }

I have skimmed the Shapefile driver and not noticed anything that should
defeat this logic but I might be missing something.  If you can demonstrate
the issue in some fashion I can easily reproduce, I'd like to dig into it.
Perhaps using:

ogrinfo -ro -spat -180 -90 180 90 ..

on your shapefile.  If you can do that, and can make the shapefile available
to me to reproduce the problem with I should be good to go.  You should
also file a Trac ticket at that point though you don't need to attach
the dataset if it is large or private.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org