text: #!/usr/bin/env python3
  ########################
  #Author: Heresh Fattahi
  #Copyright 2016
  ######################
  import argparse
  import isce
  import isceobj
  import os
  import gdal
  import numpy as np
  import xml.etree.ElementTree as ET
  
  def createParser():
      '''
      Create command line parser.
      '''
  
      parser = argparse.ArgumentParser( description='Create DEM simulation for merged images')
      parser.add_argument('-l','--lat', dest='latFile', type=str, required=True,
              help = 'latitude file in radar coordinate')
      parser.add_argument('-L','--lon', dest='lonFile', type=str, required=True,
              help = 'longitude file in radar coordinate')
      parser.add_argument('-f', '--filelist', dest='prodlist', type=str, required=True,
              help='Input file to be geocoded')
      parser.add_argument('-o', '--xoff', dest='xOff', type=int, default=0,
              help='Offset from the begining of geometry files in x direction. Default 0.0')
      parser.add_argument('-p', '--yoff', dest='yOff', type=int, default=0,
              help='Offset from the begining of geometry files in y direction. Default 0.0')
      parser.add_argument('-r', '--resampling_method', dest='resamplingMethod', type=str, default='near',
              help='Resampling method (gdalwarp resamplin methods)')
  
      
      parser.add_argument('-b', '--bbox', dest='bbox', type=str, default='',
              help='Bounding box (SNWE)')
      parser.add_argument('-x', '--lon_step', dest='lonStep', type=str, default=0.001,
              help='output pixel size (longitude) in degrees. Default 0.001')
      parser.add_argument('-y', '--lat_step', dest='latStep', type=str, default=0.001,
              help='output pixel size (latitude) in degrees. Default 0.001')
  
  
      parser.add_argument('-t', '--tiff', dest='istiff', action='store_true', default=False,
              help='Create GEOTIFF instead of standard ENVI / ISCE files')
      parser.add_argument('--alex', dest='isAlexGrid', default=False, 
              action='store_true', help='Geocode to the Antarctica grid for optical offsets used by Alex Gardner')
  
  
      return parser
  
  def cmdLineParse(iargs = None):
      '''
      Command line parser.
      '''
  
      parser = createParser()
      inps =  parser.parse_args(args = iargs)
  
  
      if not inps.isAlexGrid:
          inps.bbox = [val for val in inps.bbox.split()]
          if len(inps.bbox) != 4:
              raise Exception('Bbox should contain 4 floating point values')
          inps.outproj = 'EPSG:4326'
      else:
          print('Ignoring bbox and spacing inputs. Using standard grid for Antarctica.')
          inps.lonStep = 240.0
          inps.latStep = 240.0
          inps.outproj = 'EPSG:3031'
  
      inps.prodlist = inps.prodlist.split()
      return inps
  
  def prepare_lat_lon(inps):
  
      latFile = os.path.abspath(inps.latFile)
      lonFile = os.path.abspath(inps.lonFile)
      #cmd = 'isce2gis.py vrt -i ' + latFile
      #os.system(cmd)
      #cmd = 'isce2gis.py vrt -i ' + lonFile
      #os.system(cmd)
  
     
      width, length =  getSize(latFile)
      widthFile , lengthFile = getSize(inps.prodlist[0])
  
      print("size of lat and lon files (width, length) ", width, length)
      print("size of input file to be geocoded (width, length): ", widthFile , lengthFile)
  
      xOff = inps.xOff
      yOff = inps.yOff
  
      cmd = 'gdal_translate -of VRT -srcwin ' + str(xOff) + ' ' + str(yOff) \
             +' '+ str(width - xOff) +' '+ str(length - yOff) +' -outsize ' + str(widthFile) + \
             ' '+ str(lengthFile)  + ' -a_nodata 0 ' + latFile +'.vrt' + ' tempLAT.vrt'
  
      os.system(cmd)
  
      cmd = 'gdal_translate -of VRT -srcwin ' + str(xOff) + ' ' + str(yOff) \
            +' '+ str(int(width-xOff)) +' '+ str(int(length-yOff)) +' -outsize ' + str(widthFile) +\
             ' '+ str(lengthFile)  + ' -a_nodata 0 ' + lonFile +'.vrt' + ' tempLON.vrt'
  
      os.system(cmd)
  
      return 'tempLAT.vrt', 'tempLON.vrt' 
  
      # gdal_translate -of VRT -srcwin  384 384 64889 12785 -outsize 1013 199 ../../COMBINED/GEOM_MASTER/LAT.rdr LAT_off.vrt
      
  
  def writeVRT(infile, latFile, lonFile):
  #This function is modified from isce2gis.py
  
              #cmd = 'isce2gis.py vrt -i ' + infile
              #os.system(cmd)
  
              tree = ET.parse(infile + '.vrt')
              root = tree.getroot()
  
              meta = ET.SubElement(root, 'metadata')
              meta.attrib['domain'] = "GEOLOCATION"
              meta.tail = '\n'
              meta.text = '\n    '
  
  
              rdict = { 'Y_DATASET' : os.path.relpath(latFile , os.path.dirname(infile)),
                        'X_DATASET' :  os.path.relpath(lonFile , os.path.dirname(infile)),
                        'X_BAND' : "1",
                        'Y_BAND' : "1",
                        'PIXEL_OFFSET': "0",
                        'LINE_OFFSET' : "0",
                        'LINE_STEP' : "1",
                        'PIXEL_STEP' : "1" }
  
              for key, val in rdict.items():
                  data = ET.SubElement(meta, 'mdi')
                  data.text = val
                  data.attrib['key'] = key
                  data.tail = '\n    '
  
              data.tail = '\n'
              tree.write(infile + '.vrt')
  
  
  def runGeo(inps):
  
      #for file in inps.prodlist:
         #cmd = 'isce2gis.py envi -i ' + file
         #os.system(cmd)
  
  
      if not inps.isAlexGrid:
          WSEN = str(inps.bbox[2]) + ' ' + str(inps.bbox[0]) + ' ' + str(inps.bbox[3]) + ' ' + str(inps.bbox[1])
          latFile, lonFile = prepare_lat_lon(inps)
  
          getBound(latFile,float(inps.bbox[0]),float(inps.bbox[1]),'lat')
          getBound(lonFile,float(inps.bbox[2]),float(inps.bbox[3]),'lon')
      
          for infile in inps.prodlist:
              infile = os.path.abspath(infile)
              print ('geocoding ' + infile)
              outFile = os.path.join(os.path.dirname(infile), "geo_" + os.path.basename(infile))
              #cmd = 'isce2gis.py vrt -i '+ file + ' --lon ' + lonFile + ' --lat '+ latFile
              #os.system(cmd)
              writeVRT(infile, latFile, lonFile)
  
              cmd = 'gdalwarp -of ENVI -geoloc  -te '+ WSEN + ' -tr ' + str(inps.latStep) + ' ' + str(inps.lonStep) + ' -srcnodata 0 -dstnodata 0 ' + ' -r ' + inps.resamplingMethod + ' ' + infile +'.vrt '+ outFile
              print (cmd)
              os.system(cmd)
  
              write_xml(outFile)
  
  
      else:
          from geo2ant import getGridLimits
          ylims, xlims = getGridLimits(latfile=latFile, lonfile=lonFile)
  
          WSEN = str(xlim[0]) + ' ' + str(ylim[0]) + ' ' + str(xlim[1]) + ' ' + str(ylim[1])
          if inps.istiff:
              ext = '.tif'
              outformat = 'GTiff'
          else:
              ext = '.ant'
              outformat = 'ENVI'
  
          for infile in inps.prodlist:
              print('geocoding: ' + infile)
  
              writeVRT(infile, latFile, lonFile)
  
              cmd = 'gdalwarp -of ' + outformat + ' -t_srs ' + inps.outproj + ' -geoloc -te ' + WSEN + ' -tr ' + str(inps.lonStep) + ' ' + str(inps.latStep) + ' -srcnodata 0 -dstnodata 0 -r ' + inps.resamplingMethod + ' ' + infile + '.vrt ' + infile+'ext'
              status = os.system(cmd)
              if status:
                  raise Exception('Command {0} Failed'.format(cmd))
  
              if not inps.istiff:
                  write_xml(infile+ext)
          
  
  def getSize(infile):    
  
      ds=gdal.Open(infile + ".vrt")
      b=ds.GetRasterBand(1)
      return b.XSize, b.YSize
  
  def getBound(infile,minval,maxval,latlon): #added by Minyan Zhong
      
      ds=gdal.Open(infile)
      b=ds.GetRasterBand(1)
      data=b.ReadAsArray()
  
      idx=np.where((data>=minval) & (data<=maxval))
     
      if latlon=='lat':
          print('latitide bound in cliped area:')
      else:
          print('longitude bound in cliped area:')
      print(np.min(data[idx]),np.max(data[idx]))
  
  
  def get_lat_lon(infile):
  
      ds=gdal.Open(infile)
      b=ds.GetRasterBand(1)
      width  = b.XSize
      length = b.YSize
      minLon = ds.GetGeoTransform()[0]
      deltaLon = ds.GetGeoTransform()[1]
      maxLat = ds.GetGeoTransform()[3]
      deltaLat = ds.GetGeoTransform()[5]
      minLat = maxLat + (b.YSize)*deltaLat
  
      return maxLat, deltaLat, minLon, deltaLon, width, length
  
  def write_xml(outFile): 
  
      maxLat, deltaLat, minLon, deltaLon, width, length = get_lat_lon(outFile)
  
      unwImage = isceobj.Image.createImage()
      unwImage.setFilename(outFile)
      unwImage.setWidth(width)
      unwImage.setLength(length)
      unwImage.bands = 1
      unwImage.scheme = 'BIL'
      unwImage.dataType = 'FLOAT'
      unwImage.setAccessMode('read')
      
      unwImage.coord2.coordDescription = 'Latitude'
      unwImage.coord2.coordUnits = 'degree'
      unwImage.coord2.coordStart = maxLat 
      unwImage.coord2.coordDelta = deltaLat 
      unwImage.coord1.coordDescription = 'Longitude'
      unwImage.coord1.coordUnits = 'degree'
      unwImage.coord1.coordStart = minLon 
      unwImage.coord1.coordDelta = deltaLon 
  
     # unwImage.createImage()
      unwImage.renderHdr()
      unwImage.renderVRT()
  
  def main(iargs=None):
      '''
      Main driver.
      '''
      inps = cmdLineParse(iargs)
      runGeo(inps)
   
     
  if __name__ == '__main__':
      main()
  
  
  