Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Get extent of a raster rounded to nearest degree. Usefull in shell scripts for clipping DEMs with overlap pixels
#!/usr/bin/env python
from osgeo import gdal, ogr, osr
import sys
import logging
from optparse import OptionParser
import os
gdal.UseExceptions()
def GetExtent(gt,cols,rows):
''' Return list of corner coordinates from a geotransform
@type gt: C{tuple/list}
@param gt: geotransform
@type cols: C{int}
@param cols: number of columns in the dataset
@type rows: C{int}
@param rows: number of rows in the dataset
@rtype: C{[float,...,float]}
@return: coordinates of each corner
'''
ext=[]
xarr=[0,cols]
yarr=[0,rows]
for px in xarr:
for py in yarr:
x=gt[0]+(px*gt[1])+(py*gt[2])
y=gt[3]+(px*gt[4])+(py*gt[5])
ext.append([x,y])
yarr.reverse()
return ext
def rounded_extent(source_file, options):
try:
src_ds = gdal.Open(source_file)
except RuntimeError, e:
print 'Unable to open ' + source_file
print e
sys.exit(1)
gt = src_ds.GetGeoTransform()
cols = src_ds.RasterXSize
rows = src_ds.RasterYSize
ul, ll, lr, ur = GetExtent(gt,cols,rows)
print int(round(ll[0])), int(round(ll[1])), int(round(ur[0])), int(round(ur[1]))
if __name__=='__main__':
usage = "usage: %prog "
parser = OptionParser(usage=usage,
description="")
parser.add_option("-d", "--debug", action="store_true", dest="debug")
parser.add_option("-q", "--quiet", action="store_true", dest="quiet")
(options, args) = parser.parse_args()
logging.basicConfig(level=logging.DEBUG if options.debug else
(logging.ERROR if options.quiet else logging.INFO))
for path in args:
if not os.path.exists(path):
print path + " not found"
continue
rounded_extent(path, options)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.