[Gdal-dev] How to Create a Geometry Filter

[Roger André]

Sorry if I'm stepping on Frank here, but I'm scratching my head about
this.  Probably means that I missed the point entirely.  Seems like
the key question is, "select only features in an area bounded by 4
points."  In Python, I do this like so:

#! /usr/bin/python2.4

'''Pulls all of the zipcode polys within a bounding box and returns
the Zip Code name and centroid coords.  Spatial filtering needs min_
x, min_ y, max_x, max_y values'''

import ogr

drv = ogr.GetDriverByName( 'ESRI Shapefile' )
ds = drv.Open('zt06_d00.shp')
layer = ds.GetLayer(0)
filter_layer = layer.SetSpatialFilterRect(-122.166734286713,37.89091665754511,-122.3113385951359,37.79866675446653)
feature = layer.GetNextFeature()
name_field = feature.GetFieldIndex( 'NAME' )
while feature is not None:
  name = feature.GetField( name_field )
  print name
  polygon = feature.GetGeometryRef()
  centroid = polygon.Centroid()
  print centroid
  feature = layer.GetNextFeature()

Easy, simple, and works well for me.
--

On 10/10/07, Frank Warmerdam <warmerdam at pobox.com> wrote:
> Kent Eschenberg wrote:
> > Greetings,
> >
> > This is a combination complaint/question. I'm trying to use the layer
> > classes'  geometry filter to select only features in an area bounded by
> > 4 points.
> >
> > I can't do that directly but must use SetSpatialFilter and first create
> > an instance of OGRPolygon. I do so and look for a way to set the
> > polygon's points.
> >
> > I can't do that directly but must use addRing and first create an
> > instance of OGRLinearRing. I do that and look for a way to set the
> > ring's points.
> >
> > I can't do that directly but must use importFromWkb and first create a
> > buffer of "the well known binary data".
>
> Kent,
>
> The convenience function SetSpatialFilterRect() shows the simpliest
> approach (in C++) to create a filter polygon with four points.
>
> void OGRLayer::SetSpatialFilterRect( double dfMinX, double dfMinY,
>                                       double dfMaxX, double dfMaxY )
>
> {
>      OGRLinearRing  oRing;
>      OGRPolygon oPoly;
>
>      oRing.addPoint( dfMinX, dfMinY );
>      oRing.addPoint( dfMinX, dfMaxY );
>      oRing.addPoint( dfMaxX, dfMaxY );
>      oRing.addPoint( dfMaxX, dfMinY );
>      oRing.addPoint( dfMinX, dfMinY );
>
>      oPoly.addRing( &oRing );
>
>      SetSpatialFilter( &oPoly );
> }
>
> In short, it is not normal practice to create WKB in programs, since
> that is ... not very convenient.  It is a useful interchange format
> with other systems (ie. postgis).
>
> The OGRLinearRing isn't a "proper" OGC simple features geometry and
> so has no corresponding it WKB format.  OGRLinearRings are only
> allowed to exist in the context of a polygon.
>
> I hope this helps.
>
> 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
>
> _______________________________________________
> Gdal-dev mailing list
> Gdal-dev at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/gdal-dev
>