#! /usr/bin/python

"""Takes a point and calculates what the pixel value is of an intersecting raster."""

# USAGE: r.intersect.py <raster> <points_csv> <new_field_name>
# ASSUMES DATA IS IN THIS FORMAT
# pop|name|state|fips|lon lat

import sys
import numpy
from osgeo import gdal

def testExtents(point, ul, lr):
  # TESTS TO SEE IF POINT IS WITHIN THE EXTENTS OF THE RASTER
  point_x, point_y = point
  x_min, y_max = ul
  x_max, y_min = lr
  if not x_min <= point_x <= x_max or not y_min <= point_y <= y_max:
    return 0
  else:
    return 1

def calcIndex(point, pix_x, ul_x, pix_y, ul_y):
  # CALCULATE ARRAY POSITION OF POINT IN RASTER
  point_x, point_y = point
  x_diff = point_x - ul_x
  y_diff = ul_y - point_y
  x_pixels = abs(x_diff/pix_x)
  y_pixels = abs(y_diff/pix_y)
  x_off = int(x_pixels)
  y_off = int(y_pixels)
  offset = [x_off, y_off]
  return offset

def main():
  src_file = sys.argv[1]
  cities_file = sys.argv[2]
  new_field = sys.argv[3]
  
  # GATHER BASIC INFO ABOUT RASTER
  src_ds = gdal.Open( src_file )
  numbands = src_ds.RasterCount
  if numbands > 1:
    print "ERROR: Raster must be single-band"
    sys.exit(1)
  src_band = src_ds.GetRasterBand(1)
  src_geotransform = src_ds.GetGeoTransform()
  width = src_ds.RasterXSize
  pix_x = float(src_geotransform[1])
  ul_x = float(src_geotransform[0])
  height = src_ds.RasterYSize
  pix_y = float(src_geotransform[5])
  ul_y = float(src_geotransform[3])
  lr_x = (pix_x * width) + ul_x
  lr_y = (pix_y * height) + ul_y
  ul = [ul_x, ul_y]
  lr = [lr_x, lr_y]

  # READ POINT DATA, ASSUMES DATA IS IN THIS FORMAT
  # pop|name|state|fips|lon lat
  city_data = open(cities_file).readlines()
  print "%s|%s" % (city_data[0].strip(), new_field)
  for city in city_data[1:]:
    data_elem = city.split("|")
    lon_lat = data_elem[4].split()
    lon, lat = lon_lat
    point = (float(lon), float(lat))
    if testExtents(point, ul, lr):
      offset = calcIndex(point, pix_x, ul_x, pix_y, ul_y)
      # GET RASTER VALUE AT SPECIFIC ROW/COLUMN OFFSET IN RASTER
      iX = offset[0]
      iY = offset[1]
      src_data = src_band.ReadAsArray(0,iY,src_ds.RasterXSize,1)
      col_values = src_data[0]
      raster_value = col_values[iX]
      print "%s|%f" % (city.strip(), raster_value)
    else:
      print "%s|no_data" % (city.strip())

if __name__ == "__main__":

  main()