[FWTools] Get extent from rasters

[Trent M Hare]

I would recommend looking at gdalinfo.py (by Even Rouault)  from:
http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples

specifically the section
#/* -------------------------------------------------------------------- 
*/
#/*      Report corners. */
#/* -------------------------------------------------------------------- 


Regards,
Trent



--------------------------

"Bonus" section with more information than you needed. I already wrote it 
up so here it is. If there is a better method to do this I would like to 
hear it.

If the image has a large collar of NoDATA this can lead to unpredictable 
bounds. Using a couple gdal tricks you can get to the exact footprint (and 
bounds) by converting to a shapefile or other vector format.

1.) First map all pixels non-nodata values (values > 0) to 255.
> val_repl_greater.py  -innd 0 -outnd 255 input.tif binary_image.tif
--This routine maps all values greater than 0 to 255. Here 255 can be any 
value say if you want to tag an image with a value. e.g. Image #1, image 
#34.

The routine "val_repl_greater.py" is simply the sample "val_repl.py" (from 
link above) with one line change
replace 
      scanline = 
numpy.choose(numpy.equal(scanline,inNoData),(scanline,outNoData))
with
      scanline = 
numpy.choose(numpy.greater(scanline,inNoData),(scanline,outNoData))

2.) create polygon, The pixel value set above will be pushed into a field. 
In this case Field "DN" will have 255.
a.) If there is NO Nodata set, but the background value is 0 then you can 
"mask" using itself.
> gdal_polygonize -mask binary_image.tif  binary_image.tif  -f "ESRI 
Shapefile"  footprint1.shp

b.) if Nodata is explicitly defined, then you can skip using the mask
> gdal_polygonize binary_image.tif -f "ESRI Shapefile" footprint1.shp
note: I haven't looked into how to add in the filename into a polygon 
field which would be needed before merging. But as a work around you can 
use the index trick mentioned in step 1.

tip: here I guess you could use the sample ogrinfo.py to now report the 
polygon's extent. Which should be the raster's true pixel extent without 
the NoDATA collar.

3.) to merge many images footprints (probably should be same projection). 
#copy first file into new file
> ogr2ogr -f "ESRI Shapefile" toMerged.shp footprint1.shp
# now append remaining.
> ogr2ogr -f "ESRI Shapefile" -update -append -nln toMerged toMerged.shp 
footprint2.shp
> ogr2ogr -f "ESRI Shapefile" -update -append -nln toMerged toMerged.shp 
footprint3.shp
> ogr2ogr -f "ESRI Shapefile" -update -append -nln toMerged toMerged.shp 
footprint4.shp
...
from: http://www.gdal.org/ogr/drv_shapefile.html






From:
Susana Iraiis Delgado Rodriguez <susana.delgado_s at utzmg.edu.mx>
To:
fwtools at lists.maptools.org
Date:
10/21/2010 08:12 AM
Subject:
[FWTools] Get extent from rasters
Sent by:
fwtools-bounces at lists.maptools.org



Hello list!
 
I´m working in a python module to read rasters and get information out of 
them. I'm using GDAL and some of the osgeo libraries. I already have the 
ratsers' projection, but I also need its extent (x_min, y_min, x_max, 
y_max).
I have been searching in google and don't have any results. My code is:
 
from osgeo import gdal
from osgeo.gdalconst import *
 
gdal.AllRegister()
dataset = gdal.Open(filepath, GA_ReadOnly)
if dataset:
 _______________________________________________
FWTools mailing list
FWTools at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/fwtools
http://fwtools.maptools.org/


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/fwtools/attachments/20101021/3f1479ff/attachment.htm