14

How do I get the extent out of the GeoTiff format? I want something like this:

Extent(293518.1886150768,5680494.194041155,890338.5054657329,6267530.571271311)

AKA the xmin, xmax, ymin, ymax coordinates from the GeoTiff?

Edit: I'm implementing my own reader so I will need to know how the file format tells me this information.

Improve this question
4
  • 1
    Do you have a particular software or API in mind?
    – artwork21
    Jul 7, 2014 at 17:52
  • Sorry, I was unclear, editing the question.
    – Johan S
    Jul 7, 2014 at 17:55
  • Listgeo gives you a rather low level access to geotiff tags remotesensing.org/geotiff/listgeo.html
    – user30184
    Jul 8, 2014 at 7:24
  • hi, yes the problem is that I don't know which tags give me the information. Do you know?
    – Johan S
    Jul 8, 2014 at 9:38

5 Answers 5

Reset to default
19

Building on what @David mentioned you may use open source gdal library using python module to get image extent like this:

import gdal
from gdalconst import GA_ReadOnly

data = gdal.Open('C:/Temp/myimage.tif', GA_ReadOnly)
geoTransform = data.GetGeoTransform()
minx = geoTransform[0]
maxy = geoTransform[3]
maxx = minx + geoTransform[1] * data.RasterXSize
miny = maxy + geoTransform[5] * data.RasterYSize
print [minx, miny, maxx, maxy]
data = None

[-86.20782844487971, 40.7246942570317, -84.5446284448797, 41.8290942570317]

Reference: Find Extents of GDAL Raster

Improve this answer
1
  • Hi, thanks. As I have mentioned before I need the Keys/Tags that correspond to these values. The GeoTiff has a ModelTransformationTag, is that the transformation tag your using?
    – Johan S
    Jul 8, 2014 at 9:41
17

You can use Rasterio to get the bounding box as follows:

import rasterio

dataset = rasterio.open('example.tif')

>>> dataset.bounds
BoundingBox(left=358485.0, bottom=4028985.0, right=590415.0, top=4265115.0)
Improve this answer
0
5
+50

The tags you're interested in are: ModelTiepointTag, ModelPixelScaleTag, and ModelTransformationTag. The specification describes how they stored the information:

http://docs.opengeospatial.org/is/19-008r4/19-008r4.html#_raster_to_model_coordinate_transformation_requirements

For most common applications, the transformation between raster space and model space may be defined with a set of raster-to-model tiepoints and scaling parameters. The ModelTiepointTag and ModelPixelScaleTag may be used for this purpose.

Alternatively, the ModelTransformationTag may be used to specify the transformation matrix between the raster space (and its dependent pixel-value space) and the (possibly 3D) model space.

The ModelTiepointTag SHALL have 6 values for each of the tiepoints

The ModelPixelScaleTag SHALL have 3 values representing the scale factor in the X, Y, and Z directions

The ModelTransformationTag SHALL have 16 values representing the terms of the 4 by 4 transformation matrix. The terms SHALL be in row-major order

You could have a look at how GDAL implements them in this file:

http://trac.osgeo.org/gdal/browser/trunk/gdal/frmts/gtiff/geotiff.cpp

Improve this answer
1
  • 1
    I'm done with the implementation and these are the tags I used. What happened was that my reader didn't have to do the actual projections so the extent was all I needed, and to get them I used the tags you mentioned above.
    – Johan S
    Aug 3, 2014 at 20:07
1

Using the arcpy site package in Python, you can accomplish this by converting your geotiff to a raster object and using the extent and *max&*min classes.

import arcpy

arcpy.env.workspace = r'C:\path\to\ws'

rasters = arcpy.ListRasters()

for ras in rasters:
    f = arcpy.Raster(ras)
    xmin = f.extent.XMin
    ymin = f.extent.YMin
    xmax = f.extent.XMax
    ymax = f.extent.YMax
    rectangle = "%s %s %s %s" % (xmin, xmax, ymin, ymax)
Improve this answer
5
  • 1
    It's going to be faster if you set f equal to the extent, instead of accessing it each time. I had code similar to this in a tight loop (building raster indices), and it cut execution time in half by making that simple change. It's actually a two liner with getattr().
    – Paul
    Jul 31, 2014 at 19:17
  • Thanks for the comment @Paul. My benchmark tests indicate 0.05 seconds to extract the extent using this method. Two-liners are not necessarily better than larger sections of code--after all "Readability counts." (from the Zen of Python: import this)
    – Aaron ♦
    Jul 31, 2014 at 19:27
  • Huh, I noticed a substantial difference in my benchmarks over 1000 rasters. But for a few, yeah there won't be any noticeable difference.
    – Paul
    Jul 31, 2014 at 19:30
  • 2
    The source for import this is amusing indeed. ;) I find [getattr(extent, x) for x in ("XMin", "YMin", "XMax", "YMax")] pretty readable, but that's just me.
    – Paul
    Jul 31, 2014 at 19:36
  • 1
    I'm sorry but I want to know how the arcpy library does this, I'm writing my own reader in Scala.
    – Johan S
    Jul 31, 2014 at 19:43
0

I see that no one has mentioned yet to use GDAL's GetEnvelope(). From the Docs:

Computes and returns the bounding envelope for this geometry in the passed psEnvelope structure. This method is the same as the C function OGR_G_GetEnvelope().

So what you want to do is:

(xMin, xMax, yMin, yMax) = geom.GetEnvelope()

I feel this is more direct to what you need than using GeoTransform and doing some math to get the corners.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.