 # histogram.py

Apr 2nd, 2015
353
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import sys, struct
2. from osgeo import gdal
3. from osgeo import gdalconst
4.
5. minLat = -48
6. maxLat = -33
7. minLong = 165
8. maxLong = 179
9.
10. dataset = gdal.Open("/home/zeito/Desktop/l10g/l10g")
11. band = dataset.GetRasterBand(1)
12.
13. t = dataset.GetGeoTransform()
14. success,tInverse = gdal.InvGeoTransform(t)
15.
16. if not success:
17.     print("Failed!")
18.     sys.exit(1)
19.
20. x1, y1= gdal.ApplyGeoTransform(tInverse, minLong, minLat)
21. x2, y2= gdal.ApplyGeoTransform(tInverse, maxLong, maxLat)
22.
23. minX = int(min(x1, x2))
24. maxX = int(max(x1, x2))
25. minY = int(min(y1, y2))
26. maxY = int(max(y1, y2))
27.
28. width = (maxX - minX) + 1
29. fmt = "<" + ("h"* width)
30.
31. histogram = {}
32.
33. for y in range(minY, maxY+1):
34.     scanline = band.ReadRaster(minX, y, width, 1, width, 1,
35.                                gdalconst.GDT_Int16)
36.     values = struct.unpack(fmt, scanline)
37.
38.     for value in values:
39.         try:
40.             histogram[value] += 1
41.         except KeyError:
42.             histogram[value] = 1
43.
44. for height in sorted(histogram.keys()):
45.     print (height, histogram[height])