Advertisement
xunilk

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])
Advertisement
RAW Paste Data Copied
Advertisement