# GDAL Python cut geotiff image

I was stuck for 2 or 3 days. I have read many lectures from RS/GIS Laboratory of Utah State University. Also I read «Python Geospatial Development» book and many docs from gdal.org.

So, my problem is:

I have GeoTiff files (*.tif) with WGS84 coordinate system, (for example 35 N 532402 4892945) and need to cut this image (tif file) by specified coordinates (minX=27.37 maxX=27.42 minY=44.15 maxY=44.20) without gdalwarp utility.

So anyone have some ideas or advice?

p.s. my code for this moment:

#! /usr/bin/python
# coding: utf-8

from osgeo import gdal
from osgeo.gdalconst import *

filename = raw_input("Enter file name: ")

# Get the Imagine driver and register it
driver = gdal.GetDriverByName('GTiff')
driver.Register()

if dataset is None:
print 'Could not open ' + filename
sys.exit(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount

transform = dataset.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

for i in range(3):
x = xValues[i]
y = yValues[i]

xOffset = int((x - xOrigin) / pixelWidth)
yOffset = int((y - yOrigin) / pixelHeight)

s = str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '

for j in range(bands):
band = dataset.GetRasterBand(j+1) # 1-based index
data = band.ReadAsArray(xOffset, yOffset, 1, 1)
value = data[0,0]
s = s + str(value) + ' '
print s

• You mean you would like to do it with Python? How does your code look at the moment? – user30184 Jun 22 '16 at 13:23
• is it ok for you to use gdal_translate or is it no gdal at all? – radouxju Jun 22 '16 at 13:47
• @user30184, yep, I need to do it with Python) At this moment my code is not so pretty good) I can only read raster data with GDAL. Something like this dataset.GetGeoTransform() But also I need to convert coordinates and cut this image. And I just have no idea what to do next. – Sant Jun 22 '16 at 13:48
• please edit your question with your existing code if you want some coding help. It would be good to know why you don't want to use any gdal application , because this is the easy way. – radouxju Jun 22 '16 at 13:50
• @radouxju, 1) I cant use any console utilities(like gdalwarp or gdal_translate) 2) ok ill do this. – Sant Jun 22 '16 at 13:51

To cut an image (tif file) by using GDAL python library, and without gdalwarp utility, you need to find row and column raster indexes of top point [p1=(minX, maxY)] and bottom point [p2=(maxX, minY)]. The formulas, based in your code, are:

i1 = int((p1[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - p1[1] ) / pixelHeight)
i2 = int((p2[0] - xOrigin) / pixelWidth)
j2 = int((yOrigin - p2[1]) / pixelHeight)


Resulting (cut) raster has these new columns and rows:

new_cols = i2-i1+1
new_rows = j2-j1+1


and they can be used in the 'ReadAsArray' method for selecting the complete block of values with next instruction:

data = band.ReadAsArray(i1, j1, new_cols, new_rows)


This array can be saved easily as raster; but you need to calculate a new version of your transform because there are news xOrigin and yOrigin parameters.

I tried out my approach with next situation (points layers are used only as useful reference):

and this code:

from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
filename = "c:/pyqgis_data/aleatorio.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

print cols, rows

transform = dataset.GetGeoTransform()

print transform

p1 = (355217.199739, 4473171.2377)
p2 = (355911.113396, 4472582.9196)

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

print xOrigin, yOrigin

i1 = int((p1[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - p1[1] ) / pixelHeight)
i2 = int((p2[0] - xOrigin) / pixelWidth)
j2 = int((yOrigin - p2[1]) / pixelHeight)

print i1, j1
print i2, j2

new_cols = i2-i1+1
new_rows = j2-j1+1

data = band.ReadAsArray(i1, j1, new_cols, new_rows)

print data

new_x = xOrigin + i1*pixelWidth
new_y = yOrigin - j1*pixelHeight

print new_x, new_y

new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])

# Create gtif file
driver = gdal.GetDriverByName("GTiff")

output_file = "c:/pyqgis_data/cut_raster.tif"

dst_ds = driver.Create(output_file,
new_cols,
new_rows,
1,
gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( data )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(new_transform)

wkt = dataset.GetProjection()

# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset
dataset = None
dst_ds = None


After running the code at the Python Console of QGIS I got:

It works perfectly.

Editing Note: I use my code (sligthly modified for points and path rasters) with the same raster, but long/lat projected (EPSG: 4326), and it too worked as expected.

For this reason, I think that your issue is in this line:

pixelHeight = -transform[5]


Observe that I used a negative sign for pixelHeight.

• Ok, thx very much I have edit my code and it works but I get negative coordinates. I think I get wrong coordinates (lat/long) for cutting my images. So how to convert between Latitude/Longitude & UTM coordinates using gdal? – Sant Jun 23 '16 at 11:51
• My code works perfectly with long/lat projected rasters (see my editing note). I think that your issue is in this line: pixelHeight = -transform[5]. I used a negative value in my code. For this reason, I think that you didn't adapt correctly the complete code. – xunilk Jun 23 '16 at 13:34
• I think the issue here is that @xunik was wanting a positive value for pixelHeight even though, by convention, pixelHeight of GeoTIFF images is a negative value. So, the rest of @xunuk's code inverts every other vertical calculation to use negative numbers for example: j1 = int((yOrigin - p1[1] ) / pixelHeight), j2 = int((yOrigin - p2[1]) / pixelHeight), new_y = yOrigin - j1*pixelHeight. I suspect the code would be simpler, if pixelHeight had been left uncorrected as a negative value because then calculations involving being width and height would have the same form. – jonseymour Jun 19 at 0:39
• But, thanks @xunik, your code was a useful guide to use of the gdal API. – jonseymour Jun 19 at 0:39