# osgeo package contains the GDAl, OGR, and OSR libraries 
from osgeo import gdal
import osr
import os
import fnmatch
import sys
import string
import re
import numpy
import subprocess

# To begin, I set the initial directory workspace. "directory" is the initial directory to
# search under,
# WARNING : THIS MUST BE SET BEFORE THE USER BEGINS USING THIS SCRIPT

folder = r'C:\WorkOrder\WO_61_HUD_FFEA\hilshd_SKC_10m\img\n42w074'

# Using the os.walk and fnmatch functions, I construct a list of raster files
# to pass in my gdal.open() funtion. At the moment, I've set the file extension to '.flt'.
# WARNING : THE 'PATTERN' VARIABLE MUST BE SET BEFORE THE USER BEGINS USING
# THIS SCRIPT

pattern = "*.img"
rastList = []

# create a list of files that are within the 'folder' directory
for path, dirs, files in os.walk(folder):
	for filename in fnmatch.filter(files, pattern):
		rastList.append(os.path.join(path, filename))


# It appears that unlike OGR (which needs the ogr.GetDriverByName())
# in python GDAL detects what type of file you are opening. It loads and
# registers the correct driver, so I can just proceed to using my for loop

for raster in rastList:
       datafile = gdal.Open(raster)
       
       # command to extract information on the img size
       cols = datafile.RasterXSize
       rows = datafile.RasterYSize
       bands = datafile.RasterCount

       # access the projection information using GetProjection()
       # returns info in a WKT string format
       oldprojInfo = datafile.GetProjection()

       # create spatial refterence object for old projection

       spatialRef = osr.SpatialReference()

       spatialRef.ImportFromWkt(oldprojInfo)

       spatialRefProj = spatialRef.ExportToProj4()

       # transform the the old projection using the subprocess.call
       # and gdalwarp, could add multiprocessing...see
       # http://gis.stackexchange.com/questions/6669/converting-projected-geotiff-to-wgs84-with-gdal-and-python

       subprocess.call('gdalwarp %s %s -t_srs "+proj=lcc +lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"' %(raster, 'new_'+raster))

       # access the projection information using GetProjection()

       newprojInfo = datafile.GetProjection()

       # create a spatial refrence object for new projection

       spatialRefNew = osr.SpatialReference()

       spatialRefNew.ImportFromWkt(newprojInfo)

       spatialRefNewProj = spatialRefNew.ExportToProj4()

       # print information to shell screen

       print "Number of columns: " + str(cols)
       print "Number of rows: " + str(rows)
       print "Number of bands: " + str(bands)
       print "The original projection is (in WKT): " + oldprojInfo
       print "The original projection is (in Proj4): " + str(spatialRefProj)

       print "the new projection is (in WKT): " + newprojInfo

       print "the new projection is (in Proj4): " + str(spatialRefNewProj)
       

if datafile is None:
       print 'Could not open ' + raster
       sys.exit(1)