[gdal-dev] Extract Extent or Bounding Box for a subset of shapes for a shapefile
Hi All, I think this is the right list for the question below, as I cannot find a different list that matches my question better. I'm using ogr2ogr and ogrinfo (v1.4), and I want to extract the extent from the result of a subsetting where argument. I have a solution here that works but I'm wondering if I'm missing some more obvious solution. In a bash script, I do the following (the core snippet only, as I took out the shebang and option-getting/setting code:) --- rm -rf /tmp/tmpshpfile.* ogr2ogr -where "$WHERECLAUSE" /tmp/tmpshpfile.shp "$FILE" ogrinfo -ro -so -al /tmp/tmpshpfile.shp | \ grep Extent | \ tr -d '[:alpha:]:() -' | \ tr ',' ' ' rm -rf /tmp/tmpshpfile.* --- I can use it as follows: --- $ ./bbox.sh -f /path/to/shp/shp_counties.shp -w "COUNTY_NAME='Yellow Medicine'" 226206.484375 4934737.00313007.25 4979500.50 --- Where the argument to '-w' is put into the '$WHERECLAUSE' variable in the script. This works. As I wrote, though, am I missing something obvious in command-line switches to either of these programs that would simplify this so I wouldn't have to use this script? Or is there some other, better way of doing this? Also, is there a way to pipe the output of ogr2ogr on stdout directly to the stdin of ogrinfo? Thanks, Matt ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Extract Extent or Bounding Box for a subset of shapes for a shapefile
Hi Matthew, Your question is interesting. ogrinfo -where seemed to be the answer, but... I tried ogrinfo with the -where "sql where" option with both MapInfo and PostGIS layers. In your case it should be: ogrinfo -ro -so -al -where "$WHERECLAUSE" "whole_shapefile.shp" The returned extent is always the same, regardless of the presence or composition of the -where clause, so it refers to the whole layer rather than to the selected features. (I'm using GDAL 1.9 dev) Is this the expected behaviour of ogrinfo -where "..."? Or sould -where affect the returned Extent too? Thank you all Sig Il giorno ven, 22/04/2011 alle 08.22 -0500, Matthew Pettis ha scritto: > Hi All, > > > I think this is the right list for the question below, as I cannot > find a different list that matches my question better. > > > I'm using ogr2ogr and ogrinfo (v1.4), and I want to extract the extent > from the result of a subsetting where argument. I have a solution > here that works but I'm wondering if I'm missing some more obvious > solution. In a bash script, I do the following (the core snippet > only, as I took out the shebang and option-getting/setting code:) > > > --- > rm -rf /tmp/tmpshpfile.* > > > ogr2ogr -where "$WHERECLAUSE" /tmp/tmpshpfile.shp "$FILE" > > > ogrinfo -ro -so -al /tmp/tmpshpfile.shp | \ > grep Extent | \ > tr -d '[:alpha:]:() -' | \ > tr ',' ' ' > > > rm -rf /tmp/tmpshpfile.* > --- > > > I can use it as follows: > > > --- > $ ./bbox.sh -f /path/to/shp/shp_counties.shp -w "COUNTY_NAME='Yellow > Medicine'" > 226206.484375 4934737.00313007.25 4979500.50 > --- > > > Where the argument to '-w' is put into the '$WHERECLAUSE' variable in > the script. > > > This works. As I wrote, though, am I missing something obvious in > command-line switches to either of these programs that would simplify > this so I wouldn't have to use this script? Or is there some other, > better way of doing this? Also, is there a way to pipe the output of > ogr2ogr on stdout directly to the stdin of ogrinfo? > > > Thanks, > Matt > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Extract Extent or Bounding Box for a subset of shapes for a shapefile
Sig, ogrinfo doesn't force the recomputation of extent. It is considered as a costly operation. ogrinfo just reports the stored extent. Whereas in ogr2ogr, creating a new file automatically computes the extents. On Fri, Apr 22, 2011 at 7:44 PM, Luca Sigfrido Percich wrote: > > Hi Matthew, > > Your question is interesting. ogrinfo -where seemed to be the answer, > but... > > I tried ogrinfo with the -where "sql where" option with both MapInfo and > PostGIS layers. In your case it should be: > > ogrinfo -ro -so -al -where "$WHERECLAUSE" "whole_shapefile.shp" > > The returned extent is always the same, regardless of the presence or > composition of the -where clause, so it refers to the whole layer rather > than to the selected features. > (I'm using GDAL 1.9 dev) > > Is this the expected behaviour of ogrinfo -where "..."? Or sould -where > affect the returned Extent too? > > Thank you all > > Sig > > Il giorno ven, 22/04/2011 alle 08.22 -0500, Matthew Pettis ha scritto: > > Hi All, > > > > > > I think this is the right list for the question below, as I cannot > > find a different list that matches my question better. > > > > > > I'm using ogr2ogr and ogrinfo (v1.4), and I want to extract the extent > > from the result of a subsetting where argument. I have a solution > > here that works but I'm wondering if I'm missing some more obvious > > solution. In a bash script, I do the following (the core snippet > > only, as I took out the shebang and option-getting/setting code:) > > > > > > --- > > rm -rf /tmp/tmpshpfile.* > > > > > > ogr2ogr -where "$WHERECLAUSE" /tmp/tmpshpfile.shp "$FILE" > > > > > > ogrinfo -ro -so -al /tmp/tmpshpfile.shp | \ > > grep Extent | \ > > tr -d '[:alpha:]:() -' | \ > > tr ',' ' ' > > > > > > rm -rf /tmp/tmpshpfile.* > > --- > > > > > > I can use it as follows: > > > > > > --- > > $ ./bbox.sh -f /path/to/shp/shp_counties.shp -w "COUNTY_NAME='Yellow > > Medicine'" > > 226206.484375 4934737.00313007.25 4979500.50 > > --- > > > > > > Where the argument to '-w' is put into the '$WHERECLAUSE' variable in > > the script. > > > > > > This works. As I wrote, though, am I missing something obvious in > > command-line switches to either of these programs that would simplify > > this so I wouldn't have to use this script? Or is there some other, > > better way of doing this? Also, is there a way to pipe the output of > > ogr2ogr on stdout directly to the stdin of ogrinfo? > > > > > > Thanks, > > Matt > > ___ > > gdal-dev mailing list > > gdal-dev@lists.osgeo.org > > http://lists.osgeo.org/mailman/listinfo/gdal-dev > > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > -- Best regards, Chaitanya kumar CH. /tʃaɪθənjə/ /kʊmɑr/ +91-9494447584 17.2416N 80.1426E ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Extract Extent or Bounding Box for a subset of shapes for a shapefile
Thank you Chaitanya. I think that this makes great sense when dealing with a summary of the whole layer, especially because a lot of formats store the extent in the layer metadata/header. But if I need to dump all the features, or to get a dump or summary of a filtered set of features (in which case I presume ogr will scan all the features to apply the filter), checking for boundary limits shouldn't add much cost to the already expensive operation. Probably this goes beyond the intended scope of ogrinfo. And Matthew's workaround works fine. Thanks again Sig Il giorno ven, 22/04/2011 alle 20.20 +0530, Chaitanya kumar CH ha scritto: > Sig, > > ogrinfo doesn't force the recomputation of extent. It is considered as > a costly operation. ogrinfo just reports the stored extent. > Whereas in ogr2ogr, creating a new file automatically computes the > extents. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Extract Extent or Bounding Box for a subset of shapes for a shapefile
Luca,
You can also loop through the features with relatively simple Python
code to get the extents. See sample below (adapt at will):
def findPoints(geometry, results):
for i in range(geometry.GetPointCount()):
x,y,z = geometry.GetPoint(i)
if results['north'] == None or results['north'][1] < y:
results['north'] = (x,y)
if results['south'] == None or results['south'][1] > y:
results['south'] = (x,y)
if results['east'] == None or results['east'][0] < x:
results['east'] = (x,y)
if results['west'] == None or results['west'][0] > x:
results['west'] = (x,y)
for i in range(geometry.GetGeometryCount()):
findPoints(geometry.GetGeometryRef(i), results)
inSHP = 'test.shp'
shapefile = osgeo.ogr.Open(inSHP)
layer = shapefile.GetLayer(0)
for feat in range(layer.GetFeatureCount()):
feature = layer.GetFeature(feat)
geometry = feature.GetGeometryRef()
results = {'north' : None, 'south' : None, 'east' : None, 'west' :
None}
findPoints(geometry, results)
-marius
On Fri, 2011-04-22 at 17:39 +0200, Luca Sigfrido Percich wrote:
> Thank you Chaitanya.
>
> I think that this makes great sense when dealing with a summary of the
> whole layer, especially because a lot of formats store the extent in the
> layer metadata/header.
>
> But if I need to dump all the features, or to get a dump or summary of a
> filtered set of features (in which case I presume ogr will scan all the
> features to apply the filter), checking for boundary limits shouldn't
> add much cost to the already expensive operation.
>
> Probably this goes beyond the intended scope of ogrinfo. And Matthew's
> workaround works fine.
>
> Thanks again
>
> Sig
>
>
>
> Il giorno ven, 22/04/2011 alle 20.20 +0530, Chaitanya kumar CH ha
> scritto:
> > Sig,
> >
> > ogrinfo doesn't force the recomputation of extent. It is considered as
> > a costly operation. ogrinfo just reports the stored extent.
> > Whereas in ogr2ogr, creating a new file automatically computes the
> > extents.
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
