Merge pull request insarlab#1 from hfattahi/master

Update from original

README.txt

@@ -2,11 +2,12 @@ PySAR is an InSAR time-series package to produce three dimensional (space and ti
 
 Depending on your shell you may use commands such as the following examples to setup pysar:
 
+Using bash:
 export PYTHONPATH=/nethome/hfattahi/development/PySAR/pysar:${PYTHONPATH}
 export PATH="/nethome/hfattahi/development/PySAR:$PATH"
 export TSSARDIR=/nethome/timeseries/
-OR:
 
+Using csh:
 setenv PYTHONPATH "/nethome/hfattahi/development/PySAR" 
 set path = (/nethome/hfattahi/development/PySAR/pysar $path)
 setenv TSSARDIR "/nethome/timeseries/"

pysar/_readfile.py

@@ -23,6 +23,11 @@
 #  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 #  DEALINGS IN THE SOFTWARE.
 ###############################################################################
+#
+# Yunjun, Sep 2015: Add read_par_file()
+#                   Add read_gamma_float() and read_gamma_scomplex()
+# Yunjun, Oct 2015: Add box option for read_float32()
+
 import numpy as np
 
 def read_template(file):
@@ -48,22 +53,40 @@ def read_rsc_file(file):
   rsc_dict = dict(np.loadtxt(file,dtype=str))
   return rsc_dict
 
-def read_float32(file):
-  '''Reads roi_pac unw, cor, or hgt data.
-  
-  Requires the file path and returns amplitude and phase
-  Usage:
-    amplitude, phase, rscDictionary = readUnw('/Users/sbaker/Desktop/geo_070603-070721_0048_00018.unw')
+def read_par_file(file):
+  '''Read the .par file into a python dictionary structure.
+
+  '''
+  par_dict = dict(np.loadtxt(file,dtype=str,skiprows=2))
+  return par_dict
+
+#def read_float32(file):
+def read_float32(*args):
+  '''Reads roi_pac data (RMG format, interleaved line by line)
+        file: .unw, .cor, .hgt, .trans
+        box:  4-tuple defining the left, upper, right, and lower pixel coordinate.
+     Example:
+        a,p,r = read_float32('100102-100403.unw')
+        a,p,r = read_float32('100102-100403.unw',(100,1200,500,1500))
   '''
+
+  file = args[0]
   rscContents = read_rsc_file(file + '.rsc')
-  width = int(rscContents['WIDTH'])
+  width  = int(rscContents['WIDTH'])
   length = int(rscContents['FILE_LENGTH'])
-  oddindices = np.where(np.arange(length*2)&1)[0]
-  data = np.fromfile(file,np.float32,length*2*width)
-#  print np.shape(data)
-  data=data.reshape(length*2,width)
-  amplitude = np.array([data.take(oddindices-1,axis=0)]).reshape(length,width)
-  phase = np.array([data.take(oddindices,axis=0)]).reshape(length,width)
+  if   len(args)==1:     box = [0,0,width,length]
+  elif len(args)==2:     box = args[1]
+  else: print 'Error: only support 1/2 inputs.'; return
+
+  data = np.fromfile(file,np.float32,box[3]*2*width).reshape(box[3],2*width)
+  amplitude = data[box[1]:box[3],box[0]:box[2]]
+  phase     = data[box[1]:box[3],width+box[0]:width+box[2]]
+
+  #oddindices = np.where(np.arange(length*2)&1)[0]
+  #data = np.fromfile(file,np.float32,length*2*width).reshape(length*2,width)
+  #amplitude = np.array([data.take(oddindices-1,axis=0)]).reshape(length,width)
+  #phase     = np.array([data.take(oddindices,  axis=0)]).reshape(length,width)
+
   return amplitude, phase, rscContents
 
 def read_complex64(file):
@@ -74,13 +97,73 @@ def read_complex64(file):
     amp, phase, rscDictionary = readInt('/Users/sbaker/Desktop/geo_070603-070721_0048_00018.int')
   '''
   rscContents = read_rsc_file(file + '.rsc')
-  width = int(rscContents['WIDTH'])
+  width  = int(rscContents['WIDTH'])
   length = int(rscContents['FILE_LENGTH'])
   data = np.fromfile(file,np.complex64,length*2*width).reshape(length,width)
-  amplitude = np.array([np.hypot(data.real,data.imag)]).reshape(length,width)
-  phase = np.array([np.arctan2(data.imag,data.real)]).reshape(length,width)
+  amplitude = np.array([np.hypot(  data.real,data.imag)]).reshape(length,width)
+  phase     = np.array([np.arctan2(data.imag,data.real)]).reshape(length,width)
   return amplitude, phase, rscContents
 
+def read_gamma_float(file):
+  '''Read GAMMA FLOAT file (.mli)
+
+  '''
+  try:
+     parContents = read_rsc_file(file + '.rsc')
+     width  = int(parContents['WIDTH'])
+     length = int(parContents['FILE_LENGTH'])
+  except:
+     parContents = read_par_file(file + '.par')
+     width  = int(parContents['range_samples:'])
+     length = int(parContents['azimuth_lines:'])
+  data = np.fromfile(file,np.float32,length*width).reshape(length,width)
+  return data, parContents
+
+#def read_gamma_scomplex(file,box):
+def read_gamma_scomplex(*args):
+  '''Read GAMMA SCOMPLEX file (.slc)
+        file: complex data matrix (cpx_int16)
+        box: 4-tuple defining the left, upper, right, and lower pixel coordinate.
+     Example:
+        data,rsc = read_gamma_scomplex('100102.slc')
+        data,rsc = read_gamma_scomplex('100102.slc',(100,1200,500,1500))
+  '''
+
+  file = args[0]
+  try:
+     parContents = read_rsc_file(file + '.rsc')
+     width  = int(parContents['WIDTH'])
+     length = int(parContents['FILE_LENGTH'])
+  except:
+     parContents = read_par_file(file + '.par')
+     width  = int(parContents['range_samples:'])
+     length = int(parContents['azimuth_lines:'])
+
+  if   len(args)==1:     box = [0,0,width,length]
+  elif len(args)==2:     box = args[1]
+  else: print 'Error: only support 1/2 inputs.'; return
+  nlines = box[3]-box[1]
+  WIDTH  = box[2]-box[0]
+
+  data = np.fromfile(file,np.int16,box[3]*2*width).reshape(box[3],2*width)
+  data = data[box[1]:box[3],2*box[0]:2*box[2]].reshape(2*nlines*WIDTH,1)
+  id1 = range(0,2*nlines*WIDTH,2)
+  id2 = range(1,2*nlines*WIDTH,2)
+  real = data[id1].reshape(nlines,WIDTH)
+  imag = data[id2].reshape(nlines,WIDTH)
+
+  data_cpx = real + imag*1j
+  return data_cpx, parContents
+
+  #data = np.fromfile(file,np.int16,length*2*width).reshape(length*2,width)
+  #oddindices = np.where(np.arange(length*2)&1)[0]
+  #real = np.array([data.take(oddindices-1,axis=0)]).reshape(length,width)
+  #imag = np.array([data.take(oddindices,  axis=0)]).reshape(length,width)
+
+  #amplitude = np.array([np.hypot(  real,imag)]).reshape(length,width)
+  #phase     = np.array([np.arctan2(imag,real)]).reshape(length,width)
+  #return amplitude, phase, parContents
+
 def read_dem(file):
   '''Read a roipac dem file.
   
@@ -107,4 +190,21 @@ def read_GPS_USGS(file):
 
    return east,north,up,dates,YYYYMMDD
 
-
+def read_isce_xml(file):
+  from lxml import etree as ET
+  tree = ET.parse(file)
+  root = tree.getroot()
+  xmldict={}
+  for child in root.findall('property'):
+    attrib = child.attrib['name']
+    value = child.find('value').text
+    xmldict[attrib]=value
+  xmldict['WIDTH']=xmldict['width']
+  xmldict['FILE_LENGTH']=xmldict['length']
+  xmldict['WAVELENGTH'] = '0.05546576'
+  Date1=os.path.dirname(file).split('/')[-1].split('_')[0][2:]
+  Date2=os.path.dirname(file).split('/')[-1].split('_')[1][2:]
+  xmldict['DATE12'] = Date1 + '-' + Date2
+  xmldict['DATE1'] = Date1
+  xmldict['DATE2'] = Date2
+  return xmldict

pysar/_writefile.py

@@ -3,22 +3,40 @@
 # Copyright(c) 2013, Heresh Fattahi                        #
 # Author:  Heresh Fattahi                                  #
 ############################################################
+#
+# Yunjun, Sep 2015: Add write_gamma_float() and write_gamma_scomplex()
+# Yunjun, Oct 2015: Add support for write_float32(amp, phase, outname)
 
 import numpy as np
 
-def write_float32(data,outname):
+#def write_float32(data,outname):
+def write_float32(*args):
    # To write an array to a binary file with float32 precision
    # Format of the binary file is same as roi_pac unw, cor, or hgt data.
-
-   F=np.zeros([2*data.shape[0]*data.shape[1],1],np.float32)
-   nlines=data.shape[0]
-   WIDTH=data.shape[1]
+   # Exmaple:
+   #         write_float32(phase, outname)
+   #         write_float32(amp, phase, outname)
+
+   if len(args)==2:
+      amp     = args[0]
+      pha     = args[0]
+      outname = args[1]
+   elif len(args)==3:
+      amp     = args[0]
+      pha     = args[1]
+      outname = args[2]
+   else:
+      print 'Error while getting args: support 2/3 args only.'
+      return
+
+   nlines = pha.shape[0]
+   WIDTH  = pha.shape[1]
+   F=np.zeros([2*nlines*WIDTH,1],np.float32)
 
    for line in range(nlines):
+       F[(2*WIDTH)*(line) :       (2*WIDTH)*(line)+WIDTH]=np.reshape(amp[line][:],[WIDTH,1])
+       F[(2*WIDTH)*(line)+WIDTH : (2*WIDTH)*(line+1)]    =np.reshape(pha[line][:],[WIDTH,1])
 
-       F[(line)*(2*WIDTH):(line)*(2*WIDTH)+WIDTH]=np.reshape(data[line][:],[WIDTH,1])
-       F[(line)*(2*WIDTH)+WIDTH : (line+1)*(WIDTH*2)]=np.reshape(data[line][:],[WIDTH,1])
-
    F.tofile(outname)
 
 def write_complex64(data,outname):   
@@ -40,6 +58,21 @@ def write_dem(data,outname):
    data=np.array(data,dtype=np.int16)
    data.tofile(outname)
 
+def write_gamma_float(data,outname):
+   data=np.array(data,dtype=np.float32)
+   data.tofile(outname)
+
+def write_gamma_scomplex(data,outname):
+   # write gamma scomplex data, i.e. .slc file.
+   nlines = data.shape[0]
+   WIDTH  = data.shape[1]
+   id1 = range(0,2*nlines*WIDTH,2)
+   id2 = range(1,2*nlines*WIDTH,2)
+   F=np.zeros([2*nlines*WIDTH,1],np.int16)
+   F[id1]=np.reshape(data.real,(nlines*WIDTH,1))
+   F[id2]=np.reshape(data.imag,(nlines*WIDTH,1))
+   F.tofile(outname)
+
 
 
 

pysar/correct_dem.py

@@ -0,0 +1,50 @@
+#! /usr/bin/env python
+
+import os
+import sys
+from numpy import *
+import h5py
+from pysar import _readfile
+from pysar import _writefile
+
+def main(argv):
+  try:
+     dem_file = argv[1]
+     dem_error = argv[2]
+  except:
+    print '''
+    *******************************************
+       
+       Usage: correct_dem.py demFile geo_demErrorFile
+       Example:
+              correct_dem.py $DEMDIR/Socorro-30/Socorro_30.dem geo_DEM_error.h5
+              correct_dem.py $DEMDIR/Socorro-30/Socorro_30.dem geo_DEM_error.h5
+
+    *******************************************         
+    '''
+    sys.exit(1)
+
+
+
+  dem, demrsc = _readfile.read_dem(dem_file)
+  g = h5py.File(dem_error,'r')
+  dset  = g['dem'].get('dem')
+  dem_error = dset[0:dset.shape[0]]
+
+  print 'Correcting the DEM'
+  sum = dem + dem_error
+  print 'Creating the new DEM'
+  _writefile.write_dem(sum,'DEM_w_error.dem')
+
+  rsc_file = open('DEM_w_error.dem.rsc','w')
+  for k in demrsc.keys():
+      rsc_file.write(k+'	'+demrsc[k]+'\n')
+      rsc_file.close
+
+  date12_file=open('111111-222222_baseline.rsc','w')
+  date12_file.write('P_BASELINE_TOP_ODR'+'     '+ '000')
+  date12_file.close
+
+##########
+if __name__ == '__main__':
+   main(sys.argv[:])