5

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()

dataset = gdal.Open(filename, GA_ReadOnly)
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
9

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):

enter image description here

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:

enter image description here

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.

enter image description here

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

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.