#!/usr/bin/python
import gdal
import struct
from gdalconst import *
def interpolate(coords_list,steps):
	new_coords_x=[]
	new_coords_y=[]
	coords=[]
	for u in coords_list:
		coords.append(u.split(','))
		new_coords_x=[coords[0][0]]
		new_coords_y=[coords[0][1]]
	count=0
	for u in coords:
		change_x=0
		change_y=0
		if (count<len(coords)-1):
			for i in range(steps):
				change_x=change_x+((float(coords[count+1][0])-float(coords[count][0]))/steps)
				change_y=change_y+((float(coords[count+1][1])-float(coords[count][1]))/steps)
				x=float(coords[count][0])+change_x
				y=float(coords[count][1])+change_y
				# Need to round coords to defeats strange bug in the php-parser
				new_coords_x.append(x)
				new_coords_y.append(y)
			count=count+1
	count=0
	coords=[]
	while count<len(new_coords_x):
		coordsStr=str(new_coords_x[count])+','+str(new_coords_y[count])
		coords.append(coordsStr)
		count=count+1
	return coords
	
def r_what(layer,east_north_list):
	coords=[] #Defining list for coords
	dataset = gdal.Open(layer, GA_ReadOnly )
	geotransform = dataset.GetGeoTransform()
	band = dataset.GetRasterBand(1)
	for u in east_north_list:
		c=u.split(',')
		pix_x=int(round(((float(c[0])-geotransform[0])/geotransform[1]),0))
		pix_y=int(round(((geotransform[3]-float(c[1]))/geotransform[1]),0))
		scanline = band.ReadRaster(pix_x,pix_y, 1, 1, 1, 1, GDT_Float32 )
		z = struct.unpack('f' * 1, scanline)
		xyz=[float(c[0]),float(c[1]),z[0]]
		coords.append(xyz)
	return coords
	
def raster_query(layer,east_north,step):
	east_north_list=east_north.split(' ')
	if step:
		east_north_list=interpolate(east_north_list,int(step))
	r_what_out=r_what(layer,east_north_list)
	return r_what_out	
	
#east_north='-121.024210464,38.918681148 -119.601875424,37.748555815'
#print raster_query("/var/grassdata/wgs84/us/cellhd/S_20_05",east_north,1)