PlotVarioFunction

import os
import time as time # for reading in a timer
import numpy as np # maths functions (arrays etc.)
from osgeo import gdal, gdalconst # for reading in raster
from osgeo.gdalconst import * # for reading in raster
from matplotlib import pyplot as plt # for ploting
from scipy import signal # for convolution function
from scipy import ndimage # for resampling image
from matplotlib import cm # colour mapping
import matplotlib as matplotlib
import matplotlib.pyplot as plt
import copy as cp
from numpy import *

# start timing
startTime = time.time()

variance = []
lag = []

# Register driver
driver = gdal.GetDriverByName('ENVI')
driver.Register()

# Clear any previous plots
plt.clf()

# Set file location
file_name = r"/data/surface"

# open file
inds = gdal.Open(file_name, GA_ReadOnly)

if inds is None:
    print '\nPerhaps you need an ENVI .hdr file? If so, just open the binary up in ENVI and one will be created for you!'
    os._exit(1)
else:
    print "%s opened successfully" %file_name

print '~~~~~~~~~~~~~~'
print 'Get image size'
print '~~~~~~~~~~~~~~'
cols = inds.RasterXSize
rows = inds.RasterYSize
bands = inds.RasterCount

print "columns: %i" %cols
print "rows: %i" %rows
print "bands: %i" %bands

print '~~~~~~~~~~~~~~'
print 'Get georeference information'
print '~~~~~~~~~~~~~~'
geotransform = inds.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

print "origin x: %i" %originX
print "origin y: %i" %originY
print "width: %2.2f" %pixelWidth
print "height: %2.2f" %pixelHeight

print '~~~~~~~~~~~~~~'
print 'Convert image to 2D array'
print '~~~~~~~~~~~~~~'
band = inds.GetRasterBand(1)
image_array = band.ReadAsArray(0, 0, cols, rows)
image_array_name = file_name
print type(image_array)
print shape(image_array)

print '~~~~~~~~~~~~~~'
print 'Subsample 2D array'
print '~~~~~~~~~~~~~~'
image_array_subsample = image_array[4790:4910, 7000:7120]

for i in range (len(image_array_subsample)):
    for j in range (len(image_array_subsample[i])):
        pos_now_row = i
        pos_now_col = j
        value = image_array_subsample[i][j]
        print "%i,%i" %(pos_now_row,pos_now_col) # position
        print "value = %f" %(value)

        rows = len(image_array_subsample) # gets no. rows (count starts at 1)
        cols = len(image_array_subsample[0]) # gets no. cols (count starts at 1)
        print "Number of rows (size starting at 1): %d" %(rows)
        print "Number of columns (size starting at 1): %d" %(cols)

        rows2 = len(image_array_subsample)-1 # gets no. rows in terms of position (count starts at 0)
        cols2 = len(image_array_subsample[0])-1 # gets no. cols in terms of position (count starts at 0)
        print "Number of rows (in terms of index starting at 0): %d" %(rows2)
        print "Number of cols (in terms of index starting at 0): %d" %(cols2)

        for i2 in range (len(image_array_subsample)):
            for j2 in range (len(image_array_subsample[i2])):

                pos_now_row_2 = i2
                pos_now_col_2 = j2

                var_value = value - image_array_subsample[i2][j2]
                var_value_ABSOLUTE = math.fabs(var_value)
                variance.append(var_value_ABSOLUTE)

                lag_val_pre = ((pos_now_row_2 - pos_now_row)*(pos_now_row_2 - pos_now_row))+((pos_now_col_2 - pos_now_col)*(pos_now_col_2 - pos_now_col))
                lag_val = math.sqrt(lag_val_pre)
                lag.append(lag_val)

#### Loop through lists, set bins and create binned lists

for i in range (len(lag)):

    #print out list values
    variance_value = variance[i]
    lag_value = lag[i]
    print 'lag: %.2f variance: %.2f' %(lag_value, variance_value)

#define bins
start_bin_value = 0.0
max_lag = max(lag)
end_bin_value = max_lag
print 'max_lag: %.4f' %(max_lag)
bin_no = 20.0
print 'bin_no: %d' %(bin_no)
bin_size = (max_lag/bin_no)
print 'bin_size: %.4f' %(bin_size)

# Returns evenly spaced values within a given interval
# numpy.digitize() will "bin" the values; i.e. x_binned[i] will give the bin number of value in i-th index
# converts variance list into a numpy array
bins = np.arange(start_bin_value, end_bin_value, bin_size)
x_binned = np.digitize(lag,bins)#,right=False) <-- deprecated from version on UNIX at Swansea?
y_numpyarray = np.array(variance)

# Calc. mean var according to bin
y_means = np.array([
    (y_numpyarray[x_binned == 1].mean()),
    (y_numpyarray[x_binned == 2].mean()),
    (y_numpyarray[x_binned == 3].mean()),
    (y_numpyarray[x_binned == 4].mean()),
    (y_numpyarray[x_binned == 5].mean()),
    (y_numpyarray[x_binned == 6].mean()),
    (y_numpyarray[x_binned == 7].mean()),
    (y_numpyarray[x_binned == 8].mean()),
    (y_numpyarray[x_binned == 9].mean()),
    (y_numpyarray[x_binned == 10].mean()),
    (y_numpyarray[x_binned == 11].mean()),
    (y_numpyarray[x_binned == 12].mean()),
    (y_numpyarray[x_binned == 13].mean()),
    (y_numpyarray[x_binned == 14].mean()),
    (y_numpyarray[x_binned == 15].mean()),
    (y_numpyarray[x_binned == 16].mean()),
    (y_numpyarray[x_binned == 17].mean()),
    (y_numpyarray[x_binned == 18].mean()),
    (y_numpyarray[x_binned == 19].mean()),
    (y_numpyarray[x_binned == 20].mean())])

variance_mean_value = [(y_means[i]) for i in range(len(y_means))]
print 'number of mean values: %i' %(len(variance_mean_value))

## Loop through lag and bin arrays to cross check
for i in range(len(lag)):
    lag_value = lag[i]
    bin_value = x_binned[i]
    var_value = variance[i]
    print 'lag: %.2f, bin: %.2f, variance: %.2f' %(lag_value, bin_value, var_value)

#### plot vario-function

# make both arrays are of same type
variance_mean_np = np.array(variance_mean_value)
#variance_mean_np = np.array(y_means)
bin_np = np.array(bins)
#bin_mid_point_np = np.array(bin_mid_point)
plt.clf()
plt.plot(bin_np,variance_mean_np)
plt.xlabel('Lag (binned)')
plt.ylabel('Variance')
test_vario_plot =  r'~/Desktop/'+ 'test_vario_plot.png'
plt.savefig(test_vario_plot)