6

I am looking for an easy way to extract a raster's extent in python.

Currently, I am using gdalinfo and then parsing out the corner coordinates. The problem with this approach is that gdalinfo is not consistent in how it reports the corner coordinates. I am working with a variety of file types (hdf, tif, img, asc, etc). I need the extent to be in decimal degrees (or DMS - that I can convert to decimal degrees).

This is being used with a PostgreSQL database so I have access to any PostGIS function, gdal, or ArcGIS, although I would prefer a function that does not involve ArcGIS unless it is fast (I am trying to index a large number of files, frequently). It can be run on either Windows or Unix.

Thank you, Kevin

6

Probably the python port of gdalinfo would help you. You can see at the top of the file that all the coordinates are reported using the GDALInfoReportCorner method:

#/* -------------------------------------------------------------------- */
#/*      Report corners.                                                 */
#/* -------------------------------------------------------------------- */
print( "Corner Coordinates:" )
GDALInfoReportCorner( hDataset, hTransform, "Upper Left", \
                      0.0, 0.0 );
GDALInfoReportCorner( hDataset, hTransform, "Lower Left", \
                      0.0, hDataset.RasterYSize);
GDALInfoReportCorner( hDataset, hTransform, "Upper Right", \
                      hDataset.RasterXSize, 0.0 );
GDALInfoReportCorner( hDataset, hTransform, "Lower Right", \
                      hDataset.RasterXSize, \
                      hDataset.RasterYSize );
GDALInfoReportCorner( hDataset, hTransform, "Center", \
                      hDataset.RasterXSize/2.0, \
                      hDataset.RasterYSize/2.0 );

you can change the method itself which is implemented at the bottom:

#/************************************************************************/
#/*                        GDALInfoReportCorner()                        */
#/************************************************************************/

def GDALInfoReportCorner( hDataset, hTransform, corner_name, x, y ):


    line = "%-11s " % corner_name

#/* -------------------------------------------------------------------- */
#/*      Transform the point into georeferenced coordinates.             */
#/* -------------------------------------------------------------------- */
    adfGeoTransform = hDataset.GetGeoTransform(can_return_null = True)
    if adfGeoTransform is not None:
        dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x \
            + adfGeoTransform[2] * y
        dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x \
            + adfGeoTransform[5] * y

    else:
        line = line + ("(%7.1f,%7.1f)" % (x, y ))
        print(line)
        return False

#/* -------------------------------------------------------------------- */
#/*      Report the georeferenced coordinates.                           */
#/* -------------------------------------------------------------------- */
    if abs(dfGeoX) < 181 and abs(dfGeoY) < 91:
        line = line + ( "(%12.7f,%12.7f) " % (dfGeoX, dfGeoY ))

    else:
        line = line + ( "(%12.3f,%12.3f) " % (dfGeoX, dfGeoY ))

#/* -------------------------------------------------------------------- */
#/*      Transform to latlong and report.                                */
#/* -------------------------------------------------------------------- */
    if hTransform is not None:
        pnt = hTransform.TransformPoint(dfGeoX, dfGeoY, 0)
        if pnt is not None:
            line = line + ( "(%s," % gdal.DecToDMS( pnt[0], "Long", 2 ) )
            line = line + ( "%s)" % gdal.DecToDMS( pnt[1], "Lat", 2 ) )

    print(line)

    return True

Make it print whatever you want :)

  • Thank you so much for the reply. It looks like that might give me more control over the data then I was getting from gdalinfo. What is the hTransform variable? I have done some research and can't seem to find a simple description of what it should be. Thanks! – kguay Sep 13 '12 at 14:30
  • It is a handle to a coordinate transformation. The link that I gave you to the example is self-contained. You will see it gets created when the dataset has a coordinate system defined hTransform = osr.CoordinateTransformation( hProj, hLatLong ) – Ragi Yaser Burhum Sep 13 '12 at 17:27
0

With Arcgis 9.3 I believe it is consistent just WHEN the rasters are georeferenced!. Then after creating geoprocessing object u can reach it easily with Extent object.

extentofXmin= gp.Describe(rastersource).Extent.Xmin #this is a topleft corner

The other corners respectively reached by Ymin, Xmax, Ymax. You can see the corner situation in RasterClip in Arcgis like in this img: http://i.stack.imgur.com/aVEWV.png

  • Unfortunately, this does not give me the coordinates in decimal degrees or DMS. – kguay Sep 12 '12 at 14:21
  • True.. Somehow forgot, as for me it wasn't important the units. So, sorry - no idea how to do it fast. Just looked at arcgis 10 and seems someone has the same problem forums.arcgis.com/threads/… – najuste Sep 12 '12 at 14:40
0

I believe you've almost accomplished your goal parsing gdalinfo. The reason why you have not only decimal degrees is that your data seems to have several CRS. But gdalinfo will provide you with the information of the CRS of the given layer. Only thing you need to do is to create geometry from extent, transform it to the needed CRS (using TransformTo OGR function for example) and extract coordinates once again.

Also you may consider using gdaltindex. Maybe you will need to handle rasters with different CRS somehow, but I believe this will not be a great deal to write proper script. You will get a shp-file with rasters boundaries (and raster names and paths). Transform it to needed CRS and extract corner coordinates of the features.

0

You can compute it from the gdal API in Python without using gdalinfo, but this will not solve your coordinate system issue. Here is the first step to extract the coordinates:

try:
    from osgeo import gdal
    from osgeo import osr
except:
    import gdal
    import osr
hDataset = gdal.Open( "your_image", gdal.GA_ReadOnly )
adfGeoTransform = hDataset.GetGeoTransform(can_return_null = True)
if adfGeoTransform is not None:
    dfGeoXUL = adfGeoTransform[0] -buff
    dfGeoYUL = adfGeoTransform[3] +buff
    dfGeoXLR = adfGeoTransform[0] + adfGeoTransform[1] * hDataset.RasterXSize + adfGeoTransform[2] * hDataset.RasterYSize + buff
    dfGeoYLR = adfGeoTransform[3] + adfGeoTransform[4] * hDataset.RasterXSize + adfGeoTransform[5] * hDataset.RasterYSize - buff

then you need to project these corners into WGS84 (EPSG:4326) (warning: untested code below)

source = hDatasaet.GetSpatialRef()

target = osr.SpatialReference()
target.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(source, target)

pointUL = ogr.CreateGeometryFromWkt("POINT ("+str(dfGeoXUL)+" " +str(dfGeoYUL)+")")
pointLR = ogr.CreateGeometryFromWkt("POINT ("+str(dfGeoXLR)+" " +str(dfGeoYLR)+")")
pointUL.Transform(transform)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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