GetRasterBand() method for gdal in Python

Question

I was going through a tutorial book called Python Geospatial Development. On the chapter on using working with geospatial data in python there is an example of a script meant to handle and analysis raster data for the height values. I ran the following code:

import sys, struct
from osgeo import gdal
from osgeo import gdalconst


minLat = -48
maxLat = -33
minLong = 165
maxLong = 179

dataset = gdal.Open("l10g")
band = dataset.GetRasterBand(1)

t = dataset.GetGeoTransform()
success,tInverse = gdal.InvGeoTransform(t)
if not success:
    print("Failed!")
    sys.exit(1)

x1, y1= gdal.ApplyGeoTransform(tInverse, minLong, minLat)
x2, y2= gdal.ApplyGeoTransform(tInverse, maxLong, maxLat)

minX = int(min(x1, x2))
maxX = int(max(x1, x2))
minY = int(min(y1, y2))
maxY = int(max(y1, y2))

width = (maxX - minX) + 1
fmt = "<" + ("h"* width)

for y in range(minY, maxY+1):
    scanline = band.ReadRaster(minX, y, width, 1, width, 1,
                               gdalconst.GDT_Int16)
    values = struct.unpack(fmt, scanline)

    for values in values:
        try:
            histogram[value] += 1
        except KeyError:
            histogram[value] = 1
for height in sorted(histogram.keys()):
    print (height, hsitogram[height])

but I got the following error:

Traceback (most recent call last):
  File "D:\Python\Progies\Geospatial\l10g\histogram.py", line 12, in <module>
    band = dataset.GetRasterBand(1)
AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

Answer 1

It is necessary to clarify that "l10g" is a DEM tile from the GLOBE website (Python Geospatial Development, p.120) and is a valid raster (raw elevation data). You must also download the hdr (l10g.hdr, book p.120 and explanations of Gabor Farkas)

The l10g, l10g.hdr and histogram.py files should be in the same folder D:\Python\Progies\Geospatial\l10g\

To open "l10g", you must open your shell in this folder. If not you must adjust the path of the files

dataset = gdal.Open("l10g") 

or for example

dataset = gdal.Open("D:/Python/Progies/Geospatial/l10g/l10g")  

To verify that dataset is not None (no file open -> NoneType)

print dataset.GetGeoTransform()
(90.000000000000497, 0.0083333333329999992, 0.0, -5.000002148425331e-13, 0.0, -0.0083333333329999992)
print dataset.RasterCount()
1
print dataset.GetRasterBand(1).XSize
10800
print dataset.GetRasterBand(1).YSize
6000

Now, you can execute your script without the AttributeError: 'NoneType' object has no attribute 'GetRasterBand' error

Answer 2

This code (modified to include any datatype) uses the raster loaded directly as active Layer in QGIS (see next image):

import sys, struct
from osgeo import gdal
from osgeo import gdalconst
import os

minLat = -48
maxLat = -33
minLong = 165
maxLong = 179

layer = iface.activeLayer()
provider = layer.dataProvider()

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 
            'Int32':'i', 'Float32':'f', 'Float64':'d'}

path= provider.dataSourceUri()

(raiz, filename) = os.path.split(path)

dataset = gdal.Open(path)

band = dataset.GetRasterBand(1)

t = dataset.GetGeoTransform()
success,tInverse = gdal.InvGeoTransform(t)
if not success:
    print("Failed!")
    sys.exit(1)

x1, y1= gdal.ApplyGeoTransform(tInverse, minLong, minLat)
x2, y2= gdal.ApplyGeoTransform(tInverse, maxLong, maxLat)

minX = int(min(x1, x2))
maxX = int(max(x1, x2))
minY = int(min(y1, y2))
maxY = int(max(y1, y2))

width = (maxX - minX) + 1

BandType = gdal.GetDataTypeName(band.DataType)

print BandType

histogram = {}

for y in range(minY, maxY+1):

    scanline = band.ReadRaster(minX, 
                               y, 
                               width, 
                               1, 
                               width, 
                               1,
                               band.DataType)

    values = struct.unpack(fmttypes[BandType] * width, scanline)

    for value in values:
        try:
            histogram[value] += 1
        except KeyError:
            histogram[value] = 1

for height in sorted(histogram.keys()):
    print (height, histogram[height])

The l10g and l10g.hdr files loaded in QGIS:

It works nice (see image below) when I ran the code in the Python Console:

Editing Note: I found one error in the original post (for values in values:). Now, It works nice!

Answer 3

There are three problems with your code/approach, which I would like to help to resolve in a quick and explicit manner:

1. The raster data cannot be processed by the algorithm on its own. It needs a header file to access georeference information. Header files can be downloaded from here for the example dataset, as stated on page 120.

2. You have to declare histogram prior to filling it up with data. You can do it by adding the following line before the iteration:

   histogram = {}

3. You mistyped histogram on line 43. It isn't a big deal, but it have to be fixed for the algorithm to run.

Answer 4

As you already pointed out, you tried to use wrong paths.

The following is my fault: I input the wrong outputpath (it should be E:\pytest\testforG\L10622\L1data.tif while I input E:\pytest\testforG\resize0622\resize0622.tif, which didn't exist)

I want to write a tif file with gdal as the followings:

def writeTif(data,outputpath): if 'int8' in data.dtype.name: datatype = gdal.GDT_Byte elif 'int16' in data.dtype.name: datatype = gdal.GDT_UInt16 elif 'float32' in data.dtype.name: datatype = gdal.GDT_Float32 else: datatype = gdal.GDT_Float64

if len(data.shape) == 3:
    bands,height,width = data.shape
else:
    bands,(height,width) = 1,data.shape

driver = gdal.GetDriverByName('GTiff')

Because I made the wrong outputpath here, I got an error as:

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

dataset = driver.Create(outputpath,width,height,bands,datatype)

if(dataset != None):
    dataset.SetGeoTransform([73.445, 0.04, 0.0, 53.565, 0.0, -0.04])
if bands == 1:
    dataset.GetRasterBand(1).WriteArray(data)
else:
    for i in range(bands):
        dataset.GetRasterBand(i+1).WrtieArray(data[i])

dataset.FlushCache()
del dataset
driver = None

return None

Hope this shed some more like on the problem.