Tuesday 24 November 2015

GetRasterBand() method for gdal in Python


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
band = dataset.GetRasterBand(1)
AttributeError: 'NoneType' object has no attribute 'GetRasterBand'


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...