#!/usr/bin/env python
import ogr
import gdal
import subprocess
import sys
import tempfile
import os, os.path
import shutil

def clip(inraster, inshape, fieldname):
	# check folder "output" 
	outdir = "./output"
	if os.path.exists(outdir):
		shutil.rmtree(outdir)
	os.makedirs(outdir)
	# gathering nodata info
	raster = gdal.Open(inraster)
	nodata = raster.GetRasterBand(1).GetNoDataValue()
	ds = ogr.Open(inshape)
	lyr = ds.GetLayer(0)

	lyr.ResetReading()
	ft = lyr.GetNextFeature()

	while ft:
		# if there is a field contains the output file name
		value = ft.GetFieldAsString(fieldname)
		out_name = value.replace(' ', '_')
		outraster = outdir+'/'+out_name.strip()+".tif"    
		print inraster
		print outraster

		# the following code selects a polygon
		#subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, '-crop_to_cutline', '-cwhere', "'admin'='%s'" % country_name])
		o  = subprocess.check_output(['/usr/bin/gdalwarp', inraster, outraster, '-srcnodata', str(nodata), '-dstnodata', str(nodata), '-cutline', inshape, '-crop_to_cutline', '-cwhere', "'%s'='%s'" % (fieldname, value)])
		
		print o
		ft = lyr.GetNextFeature()

	ds = None
	raster = None

def unzip(tmpDir, inzip):
	subprocess.check_call(['/usr/bin/7z', 'x', '-o%s' % tmpDir, inzip], shell=False)
	files = [os.path.join(tmpDir,f) for f in os.listdir(tmpDir) ]
	for f in files:
		if f.endswith(".shp"):
			return f


if __name__ == '__main__':
# clip.py raster_data.tif shp.zip field_name
	tmpDir = tempfile.mkdtemp()
	shp = unzip(tmpDir, sys.argv[2])
	clip(sys.argv[1], shp, sys.argv[3] )
	shutil.rmtree(tmpDir)