2

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'
  • The dataset is 'NoneType' because "l10g" is not a valid raster. It's a directory where is histogram.py. Wich is the true raster's path? – xunilk Apr 1 '15 at 13:06
  • Thank you for your response. I was able to solve the problem. I was missing the header file for the DEM image which should be in the same directory as the DEM image file. – Ngeno Apr 3 '15 at 8:21
2

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

  • Thank you for your response. I was able to solve the problem. I was missing the header file for the DEM image which should be in the same directory as the DEM image file. Almost gave me a headache. – Ngeno Apr 3 '15 at 8:23
1

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.

enter image description here

  • Both of your observations were true. After solving the header file error I came upon the other errors which were easy to identify from the error logs. The missing header file was the main problem. – Ngeno Apr 3 '15 at 8:25
1

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:

enter image description here

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

enter image description here

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

  • This is the link to the equivalent code in the original post: http://pastebin.com/dSsUL30z. Only it must be changed the path to l10g and l10g.hdr files. – xunilk Apr 2 '15 at 10:09
  • Thanks @xunilk the missing header file was the main problem giving the non type error. But your alternative code gives me good insight as a learning developer. I will be comparing above and what I have to compare functionality and ease. Certainly, new learning tools. – Ngeno Apr 3 '15 at 8:29
0

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.

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.