from osgeo import gdal
import os

# RUN FROM ROCKET OSGEO Shell

ws = 'D:/rocketdata/helmers/DHI_web/'

rasters = os.listdir(ws)
for raster in rasters:
	if raster[-4:] == '.tif':
		print(raster)
		
		# get min/max by band
		gtif = gdal.Open(raster)
		input_files = ""
		for b in range(1,4):
			# get current band
			srcband = gtif.GetRasterBand(b)

			# Get raster statistics
			stats = srcband.GetStatistics(True, True)
			# print(stats)
			bmin = stats[0]
			bmax = stats[1]

			# Print the min, max, mean, stdev based on stats index
			print("[ STATS ] =  Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % (
				stats[0], stats[1], stats[2], stats[3]))
		
			# output raster
			out_tif = ws + 'FINAL_TILED/' + raster[:-4] + '_SCALE_band' + str(b) + '.tif'
			print(out_tif)
			input_files = input_files + out_tif + " "
			
			# min = int(bmin)
			# max = int(bmax) + 1
			
			# scale all bands 0 to 255 and tile into 512x512 blocks for faster display in geoserver
			# call1 = 'gdal_translate -ot Byte -b ' + str(b) + ' -scale -a_nodata none -co "TILED=YES" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" ' + raster + ' ' + out_tif
			call1 = 'gdal_translate -ot Byte -b ' + str(b) + ' -scale ' + str(bmin) + ' ' + str(bmax)  + ' 1 255 -a_nodata 0 -co "TILED=YES" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" ' + raster + ' ' + out_tif
			print(call1)
			os.system(call1)
			
		# make RGB composite
		rgb_img = ws + 'FINAL_TILED/' + raster[:-4] + '_SCALE_RGB.tif'
		call3 = 'C:\\Python32\\python.exe "c:\\Program Files\\GDAL\\gdal_merge.py" -o ' + rgb_img + ' -separate ' + input_files
		print(call3)
		os.system(call3)
	
	# ###  next run shell script on Macleod to do contrast enhancement
			
			#			tile into 512x512 blocks for faster display in geoserver
		rgb_img2 = ws + 'FINAL_TILED/' + raster[:-4] + '_SCALE_RGB_clip05.tif'
		rgb_img3 = ws + 'FINAL_TILED/' + raster[:-4] + '_SCALE_RGB_clip05_tiled.tif'
		call1 = 'gdal_translate -co "TILED=YES" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" ' + rgb_img2 + ' ' + rgb_img3
		print(call1)
		os.system(call1)
		
		# add overviews for even more efficient display in geoserver
		call2 = 'gdaladdo -r average ' + rgb_img3 + ' 2 4 8 16 32 64 128'
		print(call2)
		os.system(call2)