ISCE_INSAR/contrib/stack/topsStack/plotBursts.py

#!/usr/bin/env python3
#
# Author: Piyush Agram
# Copyright 2016
#
# Heresh Fattahi, updated for stack processing


import numpy as np 
import os
import isce
import isceobj
import datetime
import logging
import argparse
from isceobj.Util.ImageUtil import ImageLib as IML
from isceobj.Util.decorators import use_api
import s1a_isce_utils as ut
import glob
def createParser():
    '''
    Create command line parser.
    '''

    parser = argparse.ArgumentParser( description='Generate offset field between two Sentinel swaths')
    parser.add_argument('-i', '--reference', type=str, dest='reference', required=True, help='Path to folder')
    parser.add_argument('-k', '--kml', type=str, dest='shapefile', default=None, help='Path to kml')
    parser.add_argument('-f', '--figure', type=str, dest='figure', default=None, help='Path to output PDF')

    return parser


def cmdLineParse(iargs = None):
    '''
    Command line parser.
    '''

    parser = createParser()
    inps = parser.parse_args(args=iargs)

    return inps

def main(iargs=None):
    '''
    Merge burst products to make it look like stripmap.
    Currently will merge interferogram, lat, lon, z and los.
    '''


    inps=cmdLineParse(iargs)
    from osgeo import gdal, ogr, osr
    import matplotlib
    if inps.shapefile is not None:
        matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    import matplotlib.patches as patches


    swathList = ut.getSwathList(inps.reference)

    swathColors = ['r', 'g', 'b']
    shapeColors = ['FF0000','00FF00','0000FF']

    fig = plt.figure('Burst map')
    ax = fig.add_subplot(111,aspect='equal')
    
    tmin = None
    rmin = None
    
    xmin = 1e10
    ymin = 1e10
    xmax = -1e10
    ymax = -1e10


    if inps.shapefile is not None:
        ds = ogr.GetDriverByName('KML').CreateDataSource(inps.shapefile)
        srs = osr.SpatialReference()
        srs.SetWellKnownGeogCS('WGS84')
        layer = ds.CreateLayer('bursts', srs=srs)
        field_name = ogr.FieldDefn("Name", ogr.OFTString)
        field_name.SetWidth(16)
        layer.CreateField(field_name)
        field_name = ogr.FieldDefn("OGR_STYLE", ogr.OFTString)
        layer.CreateField(field_name)


    for swath in swathList:
        ifg = ut.loadProduct(os.path.join(inps.reference , 'IW{0}.xml'.format(swath)))
        minBurst = ifg.bursts[0].burstNumber
        maxBurst = ifg.bursts[-1].burstNumber

        if tmin is None:
            tmin = ifg.bursts[0].sensingStart
            dtime = ifg.bursts[0].azimuthTimeInterval
            rmin = ifg.bursts[0].startingRange
            drange = ifg.bursts[0].rangePixelSize


        if minBurst==maxBurst:
            print('Skipping processing of swath {0}'.format(swath))
            continue
    
        if not inps.shapefile:
            for ii in range(minBurst, maxBurst + 1):
                burst = ifg.bursts[ii-minBurst]
                x0 = np.round( (burst.startingRange - rmin)/drange)
                y0 = np.round( (burst.sensingStart - tmin).total_seconds()/ dtime)
                if ii %2 == 0:
                    style = 'solid'
                else:
                    style = 'dashdot'

                ax.add_patch( patches.Rectangle(
                    (x0,y0),
                    burst.numValidSamples,
                    burst.numValidLines,
                    edgecolor=swathColors[swath-1],
                    facecolor=swathColors[swath-1],
                    alpha=0.2,
                linestyle=style))

                xmin = min(xmin, x0)
                xmax = max(xmax, x0 + burst.numValidSamples)
                ymin = min(ymin, y0)
                ymax = max(ymax, y0 + burst.numValidLines)
        else:
            for ii in range(minBurst, maxBurst+1):
                burst = ifg.bursts[ii-minBurst]
                t0 = burst.sensingStart + datetime.timedelta(seconds = burst.firstValidLine * burst.azimuthTimeInterval)
                t1 = t0 + datetime.timedelta(seconds = burst.numValidLines * burst.azimuthTimeInterval)
                r0 = burst.startingRange + burst.firstValidSample * burst.rangePixelSize
                r1 = r0 + burst.numValidSamples * burst.rangePixelSize

                earlyNear = burst.orbit.rdr2geo(t0,r0)
                earlyFar = burst.orbit.rdr2geo(t0,r1)
                lateFar = burst.orbit.rdr2geo(t1,r1)
                lateNear = burst.orbit.rdr2geo(t1,r0)

                ring = ogr.Geometry(ogr.wkbLinearRing)
                from distutils.version import StrictVersion
                if StrictVersion(gdal.__version__) >= StrictVersion("3.0"):
                    ring.AddPoint(earlyNear[0], earlyNear[1])
                    ring.AddPoint(earlyFar[0], earlyFar[1])
                    ring.AddPoint(lateFar[0], lateFar[1])
                    ring.AddPoint(lateNear[0], lateNear[1])
                    ring.AddPoint(earlyNear[0], earlyNear[1])
                else:
                    ring.AddPoint(earlyNear[1], earlyNear[0])
                    ring.AddPoint(earlyFar[1], earlyFar[0])
                    ring.AddPoint(lateFar[1], lateFar[0])
                    ring.AddPoint(lateNear[1], lateNear[0])
                    ring.AddPoint(earlyNear[1], earlyNear[0])

                feature = ogr.Feature(layer.GetLayerDefn())
                feature.SetField('Name', 'IW{0}-{1}'.format(swath, ii))
                feature.SetField('OGR_STYLE', "PEN(c:#{0},w:8px)".format(shapeColors[swath-1]))
                feature.SetGeometry(ring)
                layer.CreateFeature(feature)
                feature = None

        
    if not inps.shapefile:   
        plt.ylim([ymin, ymax])
        plt.xlim([xmin, xmax])

        if inps.figure is not None:
            plt.savefig(inps.figure, format='pdf')
        else:
            plt.show()

    else:
        ds = None
        print('Wrote KML file: ', inps.shapefile)

if __name__ == '__main__' :
    '''
    Merge products burst-by-burst.
    '''

    main()