Initial commit

batch_gdal_translate.py

+#! /usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+Batch gdal_tranlate to extract an AOI define with -projwin
+
+Auteur: Loïc Lozac'h
+"""
+
+import os, argparse
+from subprocess import Popen, PIPE
+
+def search_files(directory='.', matchnaming='S1',extension='TIF',fictype='f'):
+    images=[]
+    matchnaming = matchnaming.upper()
+    extension = extension.upper()
+    fictype = fictype.lower()
+    for dirpath, dirnames, files in os.walk(directory):
+        if fictype != 'd':
+            for name in files:
+    #            print(os.path.join(dirpath, name) + " test")
+                if  name.upper().find(matchnaming) >= 0 and name.upper().endswith(extension):
+                    
+    #                print(os.path.join(dirpath, name) + " OK")
+                    abspath = os.path.abspath(os.path.join(dirpath, name))
+                    images.append(abspath)
+        elif fictype == 'd':
+            for dirname in dirnames:
+    #            print(os.path.join(dirpath, name) + " test")
+                if dirname.upper().find(matchnaming) >= 0 :
+                    
+    #                print(os.path.join(dirpath, name) + " OK")
+                    abspath = os.path.abspath(os.path.join(dirpath, dirname))
+                    images.append(abspath)
+        else:
+            print("search_files type error")
+            exit
+            
+    return images
+
+def process_command(cmd):
+    print("Starting : "+" ".join(cmd))
+    p = Popen(cmd, stdout=PIPE)
+#    p.wait()
+    output = p.communicate()[0]
+    if p.returncode != 0: 
+        print("process failed %d : %s" % (p.returncode, output))
+    print("#################################################")
+    
+
+
+
+if __name__ == "__main__":
+    # Make parser object
+    parser = argparse.ArgumentParser(description=
+        """
+        Batch gdal_tranlate to extract an AOI define with -projwin
+        """)
+    
+    parser.add_argument('-srcdir', action='store', required=True, help='Directory containing geotiff to process')
+    parser.add_argument('-projwin', action='store', required=True, help='Window extend in coordinates AND COMA SEPARATED, ex: -projwin ulx,uly,lrx,lry')
+    parser.add_argument('-outdir', action='store', required=True, help='Output directory')
+    
+    args=parser.parse_args()
+    
+    mvfiles=[]
+    mvfiles=search_files(args.srcdir)
+    
+    projwin=args.projwin.split(",")
+    if len(projwin) != 4:
+        print("Error: winproj argument must have 4 elements")
+        exit()
+    
+    if not os.path.exists(args.outdir):
+        os.mkdir(args.outdir)
+
+    cmdbase=["gdal_translate","-eco","-winproj",projwin[0],projwin[1],projwin[2],projwin[3]]
+    
+    for file in mvfiles :
+        destfile = os.path.join(args.outdir,os.path.basename(file)[:-4]+"_Subset.tif")
+        cmd=[]
+        cmd = [c for c in cmdbase]
+        cmd.append(file)
+        cmd.append(destfile)
+        process_command(cmd)
+        
+            
+        
+        
+            
+    
+    
+
+
+
+

organize_S1_by_satorbslice.py

+#! /usr/bin/python3
+# -*- coding: utf-8 -*-
+"""
+Script permettant de trier les images Sentinel-1 en fonction de leur plateforme,
+    de leur orbite relative et du slice number.
+
+Auteur: Loïc Lozac'h
+"""
+
+import xml.etree.ElementTree as ET
+import os, shutil, argparse
+
+
+def search_files(directory='.', matchnaming='S1',extension='dim',fictype='f'):
+    images=[]
+    matchnaming = matchnaming.upper()
+    extension = extension.upper()
+    fictype = fictype.lower()
+    for dirpath, dirnames, files in os.walk(directory):
+        if fictype != 'd':
+            for name in files:
+    #            print(os.path.join(dirpath, name) + " test")
+                if  name.upper().find(matchnaming) >= 0 and name.upper().endswith(extension):
+                    
+    #                print(os.path.join(dirpath, name) + " OK")
+                    abspath = os.path.abspath(os.path.join(dirpath, name))
+                    images.append(abspath)
+        elif fictype == 'd':
+            for dirname in dirnames:
+    #            print(os.path.join(dirpath, name) + " test")
+                if dirname.upper().find(matchnaming) >= 0 :
+                    
+    #                print(os.path.join(dirpath, name) + " OK")
+                    abspath = os.path.abspath(os.path.join(dirpath, dirname))
+                    images.append(abspath)
+        else:
+            print("search_files type error")
+            exit
+            
+    return images
+
+
+
+if __name__ == "__main__":
+    # Make parser object
+    parser = argparse.ArgumentParser(description=
+        """
+        Organize S1 calibrated files by satellite, relative orbit and slice for multi-temporal speckle filtering
+        """)
+    
+    parser.add_argument('-s1dir', action='store', required=True, help='Directory containing Calibrated Sentinel1 in DIMAP-BEAM format')
+    
+    args=parser.parse_args()
+    
+    s1dim=[]
+    s1dim=search_files(args.s1dir)
+    metas1=[]
+#     dirsel=[]
+    for file in s1dim :
+        #---metadata parsing
+        tree = ET.parse(file)
+        root = tree.getroot()
+                
+        #---metadata profile
+        balises = root.findall('Dataset_Sources/MDElem/MDElem/MDATTR')
+        metab={}
+        for balise in balises :
+            if balise.attrib['name'] == 'PRODUCT':
+                metab.setdefault('PRODUCT',balise.text)
+            elif balise.attrib['name'] == 'MISSION':
+                metab.setdefault('MISSION',"S"+balise.text[-2:])
+            elif balise.attrib['name'] == 'PASS':
+                metab.setdefault('PASS',balise.text[:3])
+            elif balise.attrib['name'] == 'REL_ORBIT':
+                metab.setdefault('REL_ORBIT',balise.text)
+            elif balise.attrib['name'] == 'slice_num':
+                metab.setdefault('slice_num',balise.text)
+        
+        outputdir = os.path.join(args.s1dir,metab['MISSION']+"_"+"ORB"+metab['REL_ORBIT']+"_"+"SL"+metab['slice_num'])
+        if not os.path.exists(outputdir):
+            os.mkdir(outputdir)
+
+        #move
+        print("copy : "+file[:-3]+"*")
+        print("to   : "+outputdir+"\n")
+        
+        shutil.move(file,outputdir)
+        shutil.move(file[:-3]+"data",outputdir)
+    
+            
+    
+    
+
+
+
+

unstack_MTSFSentinel1.py

+#! /usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+Unstack and converts to geotiff MTFS processed Sentinel-1
+
+Auteur: Loïc Lozac'h
+"""
+
+import xml.etree.ElementTree as ET
+import os, argparse, shutil
+from subprocess import Popen, PIPE
+from osgeo import gdal
+
+def search_files(directory='.', resolution='S1', extension='dim', fictype='f'):
+    images=[]
+    extension = extension.lower()
+    resolution = resolution.lower()
+    for dirpath, dirnames, files in os.walk(directory):
+        if fictype == 'f':
+            for name in files:
+    #            print(os.path.join(dirpath, name) + " test")
+                if extension and name.lower().endswith(extension) and name.lower().find(resolution) >= 0 :
+                    
+    #                print(os.path.join(dirpath, name) + " OK")
+                    abspath = os.path.abspath(os.path.join(dirpath, name))
+                    images.append(abspath)
+        elif fictype == 'd':
+            for dirname in dirnames:
+    #            print(os.path.join(dirpath, name) + " test")
+                if extension and dirname.lower().endswith(extension) and dirname.lower().find(resolution) >= 0 :
+                    
+    #                print(os.path.join(dirpath, name) + " OK")
+                    abspath = os.path.abspath(os.path.join(dirpath, dirname))
+                    images.append(abspath)
+        else:
+            print("search_files type error")
+            exit
+            
+    return images
+
+
+def process_command(cmd):
+    print("Starting : "+" ".join(cmd))
+    p = Popen(cmd, stdout=PIPE)
+#    p.wait()
+    output = p.communicate()[0]
+    if p.returncode != 0: 
+        print("process failed %d : %s" % (p.returncode, output))
+    print("#################################################")
+
+if __name__ == "__main__":
+    # Make parser object
+    parser = argparse.ArgumentParser(description=
+        """
+        Unstack, converts to geotiff and mosaic MTFS processed Sentinel-1 BEAM-DIMAP stack files
+        """)
+    
+    parser.add_argument('-dimdir', action='store', required=True, help='Directory containing BEAM-DIMAP stack files to be processed (.dim)')
+    parser.add_argument('-outdir', action='store', required=True, help='Output directory')
+    
+    args=parser.parse_args()
+    
+    s1dims=search_files(args.dimdir)
+    
+    if len(s1dims)==0:
+        print("No dim file found in directory!")
+        exit()
+    
+    cmdbase = ["gdal_translate"] 
+    
+    if not os.path.exists(args.outdir):
+        os.mkdir(args.outdir)
+    
+    histlist=[]
+    tomosa=[]
+    i=0
+    for s1dim in s1dims:
+        
+        tree = ET.parse(s1dim)
+        root = tree.getroot()
+
+        datadir = s1dim[:-3]+"data"
+
+        abstractnode = root.find(".//MDElem[@name='Abstracted_Metadata']")
+        mstproduct = abstractnode.find(".//*[@name='PRODUCT']").text
+        mstproductsp = mstproduct.split("_")
+    
+        slavenode = root.find(".//MDElem[@name='Slave_Metadata']")
+        
+        
+        for child in slavenode:
+            
+            if child.get("name")=="Master_bands":
+                mstbands = child.text.split(" ")
+                for b in mstbands:
+                    cmd = []
+                    i+=1
+                    bsp = b.split("_")
+                    destfile = os.path.join(args.outdir,"_".join(mstproductsp[:6])+"_"+"_".join(bsp[:2])+".tif")
+                    sourcefile = os.path.join(datadir,b+"_db.img")
+                    cmd = ["gdal_translate"]
+                    cmd.append(sourcefile)
+                    cmd.append(destfile)
+                    if not os.path.exists(destfile):
+                        process_command(cmd)
+                    histlist.append(destfile)
+                    
+                    acqdate = mstproductsp[4][:9]
+                    polari = bsp[1]
+                    plateform = mstproductsp[0]
+                    l = [histlist.index(i) for i in histlist if (plateform in i) and (acqdate in i) and (polari in i)]
+                    if len(l) == 2 :
+                        tomosa.append((histlist[l[0]],histlist[l[1]]))
+            
+            else:
+                product = child.find(".//*[@name='PRODUCT']").text
+                productsp = product.split("_")
+                bands = child.find(".//*[@name='Slave_bands']").text.split(" ")
+                for b in bands:
+                    cmd = []
+                    i+=1
+                    bsp = b.split("_")
+                    destfile = os.path.join(args.outdir,"_".join(productsp[:6])+"_"+"_".join(bsp[:2])+".tif")
+                    sourcefile = os.path.join(datadir,b+"_db.img")
+                    cmd = ["gdal_translate"]
+                    cmd.append(sourcefile)
+                    cmd.append(destfile)
+                    
+                    if not os.path.exists(destfile):
+                        process_command(cmd)
+                    
+                    histlist.append(destfile)
+                    acqdate = productsp[4][:9]
+                    polari = bsp[1]
+                    plateform = productsp[0]
+                    l = [histlist.index(i) for i in histlist if (plateform in i) and (acqdate in i) and (polari in i)]
+                    if len(l) == 2 :
+                        tomosa.append((histlist[l[0]],histlist[l[1]]))
+                    
+            
+            print("\nNumber of bands processed: "+str(i)+"\n")
+        
+    print("Mosaicing...\n")
+    
+    tempdir = os.path.join(args.outdir,"temp")
+    if not os.path.exists(tempdir):
+        os.mkdir(tempdir)
+        
+    for mos in tomosa:
+        
+        tempmos0 = os.path.join(tempdir,os.path.basename(mos[0]))
+        tempmos1 = os.path.join(tempdir,os.path.basename(mos[1]))
+        shutil.move(mos[0], tempmos0)
+        shutil.move(mos[1], tempmos1)
+        outmos = mos[0]
+        cmd = ["gdal_merge.py","-n","0","-a_nodata","0","-o",outmos,tempmos0,tempmos1]
+        print("\nOutfile: "+outmos)
+        process_command(cmd)
+            
+            
+    
+        
+            
+    
+    
+
+
+