[Gdal-dev] Re: OGR C API question

[Etienne Dube]

Hi Stephan,

You have to make sure the SRID reference from the layer table is present in 
the spatial_ref_sys table. Otherwise it seems that OGR is unable to apply the 
spatial filter. First, if you haven't actually installed the PostGIS 
functions in your database, do it:

createlang -U postgres plpgsql <dbname>
psql postgres -d <dbname> -f postgis.sql

Then, you can use ogr2ogr to fill in the rows from a shapefile:

ogr2ogr -f "PostgreSQL" -s_srs 
'GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]]'  
"PG:dbname=<dbname> user=<username>" SHAPEFILE.SHP

(Of course, you have to replace <dbname> and <username> by the correct 
values.)

I've found that if you don't specify the -s_srs (source spatial ref. sys.) 
when importing an ESRI shapefile with ogr2ogr, the SRID row in 
spatial_ref_sys won't be created. Make sure you use the same datum as the 
data in your shapefile.

Now, spatial filtering should work with OGR on this dataset. However, (correct 
me if I'm wrong) the current implementation seems to be broken, it won't use 
spatial queries to fetch the features from the table. Instead, it SELECTs all 
the data in the layer, and uses a brute force approach (intersection on all 
geometries). See the "OGRFeature *OGRPGLayer::GetNextFeature()" function in 
ogrpglayer.cpp. So, if you're working on large datasets, methods like 
OGRLayer::GetExtent() or OGRLayer::GetFeatureCount() will take a looooooong 
time to execute. Same thing if, say, you want to display all the geometries 
in a small rectangle; OGR will have to intersect this rectangle with all the 
geometries in the layer, which takes some time to do.

Etienne

On Monday 29 March 2004 10:44, Stephan Holl wrote:
> On Mon, 22 Mar 2004 23:38:46 +0100
>
> Markus Neteler <neteler at itc.it> wrote:
> > On Mon, Mar 22, 2004 at 09:53:21AM -0500, Frank Warmerdam wrote:
> > > Markus Neteler wrote:
> > > >	poSpatialFilter = OGR_G_CreateGeometry( wkbPolygon );
> > > >	Ogr_oRing = OGR_G_CreateGeometry( wkbLineString );
> > > >        OGR_G_AddPoint(Ogr_oRing, xmin, ymin, 0);
> > > >        OGR_G_AddPoint(Ogr_oRing, xmin, ymax, 0);
> > > >        OGR_G_AddPoint(Ogr_oRing, xmax, ymax, 0);
> > > >        OGR_G_AddPoint(Ogr_oRing, xmax, ymin, 0);
> > > >        OGR_G_AddPoint(Ogr_oRing, xmin, ymin, 0);
> > > >        OGR_G_AddGeometryDirectly(poSpatialFilter, Ogr_oRing);
> > > >
> > > >        OGR_L_SetSpatialFilter(Ogr_layer, poSpatialFilter );
> > > >
> > > >Does anyone see why this doesn't work any longer? The spatial
> > > >filter is simply ignored.
> > >
> > > Markus,
> > >
> > > The wkbLineString should be wkbLinearRing. It looks like the
> > > AddGeometryDirectly() method returns
> > > OGRERR_UNSUPPORTED_GEOMETRY_TYPE if you pass anything other than a
> > > wkbLinearRing to add to a polygon.   The checking may have been
> > > added after the initial implementation though I see nothing obvious
> > > in the change logs.
> >
> > Great, now it works smoothly (tried with SHAPE, hope it does as well
> > with PostGIS data).
>
> After some testin it does not work with PostGIS-data. It always pulls
> the whole geometries out of PostGIS.
>
> Does anybody know a solution?!
>
> Thanks