KJ135
/
ISCE_INSAR
mirror of http://172.16.0.12:5000/WBJY/ISCE_INSAR.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #103 from yunjunz/mask 

stripmapStack/MaskAndFitler: more plotting options
LT1AB
piyushrpt 2020-03-23 21:47:54 -07:00 committed by GitHub
parent 3f82b2cfe9 17b22417fd
commit e0b7ad05a1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 160 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

281
contrib/stack/stripmapStack/MaskAndFilter.py
View File

@ -4,52 +4,80 @@
# Copyright 2016 # Copyright 2016
# #
import numpy as np
import argparse
import os import os
import isce import argparse
import isceobj import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
import gdal import gdal
from gdalconst import GA_ReadOnly from gdalconst import GA_ReadOnly
from scipy import ndimage
# suppress the DEBUG message
import logging
mpl_logger = logging.getLogger('matplotlib')
mpl_logger.setLevel(logging.WARNING)
import isce
import isceobj
from isceobj.Util.ImageUtil import ImageLib as IML
GDAL2NUMPY_DATATYPE = { GDAL2NUMPY_DATATYPE = {
1 : np.uint8,
1 : np.uint8, 2 : np.uint16,
2 : np.uint16, 3 : np.int16,
3 : np.int16, 4 : np.uint32,
4 : np.uint32, 5 : np.int32,
5 : np.int32, 6 : np.float32,
6 : np.float32, 7 : np.float64,
7 : np.float64, 10: np.complex64,
10: np.complex64, 11: np.complex128,
11: np.complex128,
} }
EXAMPLE = '''example:
MaskAndFilter.py -d offset.bip -s offset_snr.bip
MaskAndFilter.py -d offset.bip -s offset_snr.bip --plot
'''
def createParser(): def createParser():
''' '''
Command line parser. Command line parser.
''' '''
parser = argparse.ArgumentParser( description='filters the densOffset, oversamples it and adds back to the geometry offset') parser = argparse.ArgumentParser(description='Mask and filter the densOffset',
formatter_class=argparse.RawTextHelpFormatter,
epilog=EXAMPLE)
parser.add_argument('-d', '--dense_offset', dest='denseOffset', type=str, required=True, parser.add_argument('-d', '--dense_offset', dest='denseOffset', type=str, required=True,
help='The dense offsets file obtained from cross correlation or any other approach') help='The dense offsets file obtained from cross correlation or any other approach')
parser.add_argument('-s', '--snr', dest='snr', type=str, required=True, parser.add_argument('-s', '--snr', dest='snr', type=str, required=True,
help='The SNR of the dense offsets obtained from cross correlation or any other approach') help='The SNR of the dense offsets obtained from cross correlation or any other approach')
parser.add_argument('-n', '--filter_size', dest='filterSize', type=int, default=8, parser.add_argument('-n', '--filter_size', dest='filterSize', type=int, default=8,
help='The size of the median filter') help='Size of the median filter (default: %(default)s).')
parser.add_argument('-t', '--snr_threshold', dest='snrThreshold', type=float, default=5, parser.add_argument('-t', '--snr_threshold', dest='snrThreshold', type=float, default=5,
help='The snr threshold used to mask the offset') help='Min SNR used in the offset (default: %(default)s).')
# output
parser.add_argument('-A', '--output_azimuth_offset', dest='outAzimuth', type=str, default='filtAzimuth.off', parser.add_argument('-A', '--output_azimuth_offset', dest='outAzimuth', type=str, default='filtAzimuth.off',
help='The azimuth offsets after rubber sheeting') help='File name of the azimuth offsets after rubber sheeting (default: %(default)s).')
parser.add_argument('-R', '--output_range_offset', dest='outRange', type=str, default='filtRange.off', parser.add_argument('-R', '--output_range_offset', dest='outRange', type=str, default='filtRange.off',
help='The range offsets after rubber sheeting') help='File name of the range offsets after rubber sheeting (default: %(default)s).')
parser.add_argument('-o', '--output_directory', dest='outDir', type=str, default='./', parser.add_argument('-o', '--output_directory', dest='outDir', type=str, default='./',
help='Output directory') help='Output directory (default: %(default)s).')
parser.add_argument('-p', '--plot', dest='plot', action='store_true', default=False,
help='plot the offsets before and after masking and filtering') # plot
plot = parser.add_argument_group('plot')
plot.add_argument('-p', '--plot', dest='plot', action='store_true', default=False,
help='plot the offsets before and after masking and filtering')
plot.add_argument('-v', dest='vlim', nargs=2, type=float, default=(-0.05, 0.05),
help='display range for offset (default: %(default)s).')
plot.add_argument('--v-snr', dest='vlim_snr', nargs=2, type=float, default=(0, 100),
help='display range for offset SNR (default: %(default)s).')
plot.add_argument('--figsize', dest='figsize', nargs=2, type=float, default=(18, 5),
help='figure size in inch (default: %(default)s).')
plot.add_argument('--save', dest='fig_name', type=str, default=None,
help='save figure as file')
return parser return parser
@ -58,7 +86,7 @@ def cmdLineParse(iargs = None):
return parser.parse_args(args=iargs) return parser.parse_args(args=iargs)
def read(file, processor='ISCE' , bands=None , dataType=None): def read(file, processor='ISCE', bands=None, dataType=None):
''' raeder based on GDAL. ''' raeder based on GDAL.
Args: Args:
@ -73,10 +101,13 @@ def read(file, processor='ISCE' , bands=None , dataType=None):
Returns: Returns:
* data : A numpy array with dimensions : number_of_bands * length * width * data : A numpy array with dimensions : number_of_bands * length * width
''' '''
# generate ENVI hdr file and fix the file path in xml
file = os.path.abspath(file)
if processor == 'ISCE': if processor == 'ISCE':
cmd = 'isce2gis.py envi -i ' + file img, dataname, metaname = IML.loadImage(file)
os.system(cmd) img.filename = file
img.setAccessMode('READ')
img.renderHdr()
dataset = gdal.Open(file,GA_ReadOnly) dataset = gdal.Open(file,GA_ReadOnly)
@ -97,23 +128,22 @@ def read(file, processor='ISCE' , bands=None , dataType=None):
# Fill the array with the Raster bands # Fill the array with the Raster bands
idx=0 idx=0
for i in bands: for i in bands:
band=dataset.GetRasterBand(i) band=dataset.GetRasterBand(i)
data[idx,:,:] = band.ReadAsArray() data[idx,:,:] = band.ReadAsArray()
idx+=1 idx+=1
dataset = None dataset = None
return data return data
def write(raster, fileName, nbands, bandType): def write(raster, fileName, nbands, bandType):
############
# Create the file # Create the file
driver = gdal.GetDriverByName( 'ENVI' ) driver = gdal.GetDriverByName( 'ENVI' )
dst_ds = driver.Create(fileName, raster.shape[1], raster.shape[0], nbands, bandType ) dst_ds = driver.Create(fileName, raster.shape[1], raster.shape[0], nbands, bandType )
dst_ds.GetRasterBand(1).WriteArray( raster, 0 ,0 ) dst_ds.GetRasterBand(1).WriteArray( raster, 0 ,0 )
dst_ds = None dst_ds = None
return
def fill(data, invalid=None): def fill(data, invalid=None):
""" """
@ -132,40 +162,46 @@ def fill(data, invalid=None):
if invalid is None: invalid = np.isnan(data) if invalid is None: invalid = np.isnan(data)
ind = ndimage.distance_transform_edt(invalid, ind = ndimage.distance_transform_edt(invalid,
return_distances=False, return_distances=False,
return_indices=True) return_indices=True)
return data[tuple(ind)] return data[tuple(ind)]
def mask_filter(inps, band, outName, plot=False):
#masking and Filtering
Offset = read(inps.denseOffset, bands=band)
Offset = Offset[0,:,:]
def mask_filter(inps, band, outName):
"""masking and Filtering"""
# read offset
offset = read(inps.denseOffset, bands=band)
offset = offset[0,:,:]
# read SNR
snr = read(inps.snr, bands=[1]) snr = read(inps.snr, bands=[1])
snr = snr[0,:,:] snr = snr[0,:,:]
snr[np.isnan(snr)] = 0
# Masking the dense offsets based on SNR # mask the offset based on SNR
print ('masking the dense offsets with SNR threshold: ', inps.snrThreshold) print('masking the dense offsets with SNR threshold: {}'.format(inps.snrThreshold))
Offset[snr<inps.snrThreshold] = 0 offset1 = np.array(offset)
offset1[snr < inps.snrThreshold] = np.nan
############ # fill the hole in offset with nearest data
Offset[snr<inps.snrThreshold]=np.nan print('fill the masked out region with nearest data')
Offset = fill(Offset) offset2 = fill(offset1)
############
# Median filtering the masked offsets # median filtering
print ('Filtering with median filter with size : ', inps.filterSize) print('filtering with median filter with size : ', inps.filterSize)
Offset = ndimage.median_filter(Offset, size=inps.filterSize) offset3 = ndimage.median_filter(offset2, size=inps.filterSize)
length, width = Offset.shape length, width = offset3.shape
# writing the masked and filtered offsets to a file # write data to file
print ('writing masked and filtered offsets to: ', outName) print('writing masked and filtered offsets to: ', outName)
write(Offset, outName, 1, 6) write(offset3, outName, 1, 6)
# write the xml file # write the xml/vrt/hdr file
img = isceobj.createImage() img = isceobj.createImage()
img.setFilename(outName) img.setFilename(outName)
img.setWidth(width) img.setWidth(width)
img.setLength(length)
img.setAccessMode('READ') img.setAccessMode('READ')
img.bands = 1 img.bands = 1
img.dataType = 'FLOAT' img.dataType = 'FLOAT'
@ -175,103 +211,72 @@ def mask_filter(inps, band, outName, plot=False):
img.renderVRT() img.renderVRT()
#img.finalizeImage() #img.finalizeImage()
################################ return [snr, offset, offset1, offset2, offset3]
if plot:
import matplotlib.pyplot as plt
fig = plt.figure()
ax=fig.add_subplot(1,2,1)
# cax=ax.imshow(azOffset[800:,:], vmin=-2, vmax=4)
cax=ax.imshow(Offset, vmin=-2, vmax=4)
ax.set_title('''Offset''')
ax=fig.add_subplot(1,2,2)
#ax.imshow(azOffset_filt[800:,:], vmin=-2, vmax=4)
ax.imshow(Offset_filt, vmin=-2, vmax=4)
ax.set_title('''Offset filt''')
plt.show()
def getShape(file): def plot_mask_and_filtering(az_list, rg_list, inps=None):
dataset = gdal.Open(file,GA_ReadOnly)
return dataset.RasterYSize, dataset.RasterXSize
def resampleOffset(maskedFiltOffset, geometryOffset, resampledOffset, outName): fig, axs = plt.subplots(nrows=2, ncols=5, figsize=inps.figsize, sharex=True, sharey=True)
titles = ['SNR', 'offset', 'offset (mask)', 'offset (mask/fill)', 'offset (mask/fill/filter)']
length, width = getShape(geometryOffset)
print('oversampling the filtered and masked offsets to the width and length:', width, ' ', length )
cmd = 'gdal_translate -of ENVI -outsize ' + str(width) + ' ' + str(length) + ' ' + maskedFiltOffset + ' ' + resampledOffset
os.system(cmd)
img = isceobj.createImage() # plot SNR
img.setFilename(resampledOffset) im0 = axs[0,0].imshow(az_list[0], vmin=inps.vlim_snr[0], vmax=inps.vlim_snr[1], cmap='RdBu')
img.setWidth(width) im0 = axs[1,0].imshow(rg_list[0], vmin=inps.vlim_snr[0], vmax=inps.vlim_snr[1], cmap='RdBu')
img.setAccessMode('READ') axs[0,0].set_title('SNR', fontsize=12)
img.bands = 1
img.dataType = 'FLOAT'
img.scheme = 'BIP'
#img.createImage()
img.renderHdr()
img.renderVRT()
#img.finalizeImage()
print ('Adding the dense offsets to the geometry offsets. Output: ', outName) # label
cmd = "imageMath.py -e='a+b' -o " + outName + " -t float --a=" + geometryOffset + " --b=" + resampledOffset axs[0,0].set_ylabel('azimuth', fontsize=12)
os.system(cmd) axs[1,0].set_ylabel('range', fontsize=12)
# plot offset
for i in range(1,len(az_list)):
im1 = axs[0,i].imshow(az_list[i], vmin=inps.vlim[0], vmax=inps.vlim[1], cmap='jet')
im1 = axs[1,i].imshow(rg_list[i], vmin=inps.vlim[0], vmax=inps.vlim[1], cmap='jet')
axs[0,i].set_title(titles[i], fontsize=12)
fig.tight_layout()
# colorbar
fig.subplots_adjust(bottom=0.15)
cax0 = fig.add_axes([0.09, 0.1, 0.08, 0.015])
cbar0 = plt.colorbar(im0, cax=cax0, orientation='horizontal')
cax0.yaxis.set_ticks_position('left')
#fig.subplots_adjust(right=0.93)
cax1 = fig.add_axes([0.57, 0.1, 0.08, 0.015])
cbar1 = plt.colorbar(im1, cax=cax1, orientation='horizontal')
cbar1.set_label('pixel', fontsize=12)
# save figure to file
if inps.fig_name is not None:
inps.fig_name = os.path.abspath(inps.fig_name)
print('save figure to file {}'.format(inps.fig_name))
plt.savefig(inps.fig_name, bbox_inches='tight', transparent=True, dpi=300)
plt.show()
return
def main(iargs=None): def main(iargs=None):
inps = cmdLineParse(iargs) inps = cmdLineParse(iargs)
#######################
# working on the azimuth offsets
#######################
#cmd = 'isce2gis.py envi -i ' + inps.geometryAzimuthOffset
#os.system(cmd)
if not os.path.exists(inps.outDir):
os.makedirs(inps.outDir)
####################### if not os.path.exists(inps.outDir):
# masking the dense offsets based on SNR and median filter the masked offsets os.makedirs(inps.outD)
maskedFiltOffset = os.path.join(inps.outDir, inps.outAzimuth)
mask_filter(inps, band=[1], outName = maskedFiltOffset, plot=inps.plot)
cmd = 'isce2gis.py envi -i ' + maskedFiltOffset #######################
os.system(cmd) # masking the dense offsets based on SNR and median filter the masked offs
#######################
# resampling the masked and filtered dense offsets to the same
# grid size of the geometry offsets and adding back to the
# geometry offsets
# outAzimuth = os.path.join(inps.outDir, inps.outAzimuth)
# resampledDenseOffset = os.path.join(inps.outDir, 'filtAzOff_resampled.bil') # azimuth offsets
# resampleOffset(maskedFiltOffset, inps.geometryAzimuthOffset, resampledDenseOffset, outAzimuth) inps.outAzimuth = os.path.join(inps.outDir, inps.outAzimuth)
az_list = mask_filter(inps, band=[1], outName=inps.outAzimuth)
#######################
#######################
# working on the range offsets
#######################
#cmd = 'isce2gis.py envi -i ' + inps.geometryRangeOffset
#os.system(cmd)
if not os.path.exists(inps.outDir):
os.makedirs(inps.outDir)
#######################
# masking the dense offsets based on SNR and median filter the masked offsets
maskedFiltOffset = os.path.join(inps.outDir, inps.outRange) # range offsets
mask_filter(inps, band=[2], outName = maskedFiltOffset, plot=inps.plot) inps.outRange = os.path.join(inps.outDir, inps.outRange)
rg_list = mask_filter(inps, band=[2], outName=inps.outRange)
cmd = 'isce2gis.py envi -i ' + maskedFiltOffset # plot result
os.system(cmd) if inps.plot:
####################### plot_mask_and_filtering(az_list, rg_list, inps)
# resampling the masked and filtered dense offsets to the same grid size of the geometry offsets return
#outRange = os.path.join(inps.outDir, inps.outRange)
#resampledDenseOffset = os.path.join(inps.outDir, 'filtRngOff_resampled.bil')
#resampleOffset(maskedFiltOffset, inps.geometryRangeOffset, resampledDenseOffset, outRange)
if __name__ == '__main__': if __name__ == '__main__':

13
contrib/stack/stripmapStack/createWaterMask.py
View File

@ -32,10 +32,16 @@ def createParser():
'and (0,360) or (-180,180) for longitudes.') 'and (0,360) or (-180,180) for longitudes.')
parser.add_argument('-d','--dem_file', dest='demName', type=str, default=None, parser.add_argument('-d','--dem_file', dest='demName', type=str, default=None,
help='DEM file in geo coordinates, i.e. demLat*.dem.wgs84.') help='DEM file in geo coordinates, i.e. demLat*.dem.wgs84.')
parser.add_argument('-l', '--lat_file', dest='latName', type=str, default=None, parser.add_argument('-l', '--lat_file', dest='latName', type=str, default=None,
help='pixel by pixel lat file in radar coordinate') help='pixel by pixel lat file in radar coordinate')
parser.add_argument('-L', '--lon_file', dest='lonName', type=str, default=None, parser.add_argument('-L', '--lon_file', dest='lonName', type=str, default=None,
help='pixel by pixel lat file in radar coordinate') help='pixel by pixel lat file in radar coordinate')
parser.add_argument('--fill', dest='fillValue', type=int, default=-1, choices={-1,0},
help='fill value for pixels with missing data. Default: -1.\n'
'-1 for water body\n'
' 0 for land')
parser.add_argument('-o', '--output', dest='outfile', type=str, parser.add_argument('-o', '--output', dest='outfile', type=str,
help='output filename of water mask in radar coordinates') help='output filename of water mask in radar coordinates')
return parser return parser
@ -73,7 +79,7 @@ def dem2bbox(dem_file):
return bbox return bbox
def download_waterMask(bbox, dem_file): def download_waterMask(bbox, dem_file, fill_value=-1):
out_dir = os.getcwd() out_dir = os.getcwd()
# update out_dir and/or bbox if dem_file is input # update out_dir and/or bbox if dem_file is input
if dem_file: if dem_file:
@ -86,8 +92,7 @@ def download_waterMask(bbox, dem_file):
#inps.waterBodyGeo = sw.defaultName(inps.bbox) #inps.waterBodyGeo = sw.defaultName(inps.bbox)
sw.outputFile = os.path.join(out_dir, sw.defaultName(bbox)) sw.outputFile = os.path.join(out_dir, sw.defaultName(bbox))
sw._noFilling = False sw._noFilling = False
sw._fillingValue = -1.0 #fill pixels without DEM data with value of -1, same as water body sw._fillingValue = fill_value
#sw._fillingValue = 0.0
sw.stitch(bbox[0:2], bbox[2:]) sw.stitch(bbox[0:2], bbox[2:])
return sw.outputFile return sw.outputFile
@ -106,7 +111,7 @@ def geo2radar(geo_file, rdr_file, lat_file, lon_file):
def main(iargs=None): def main(iargs=None):
inps = cmdLineParse(iargs) inps = cmdLineParse(iargs)
geo_file = download_waterMask(inps.bbox, inps.demName) geo_file = download_waterMask(inps.bbox, inps.demName, inps.fillValue)
if inps.latName and inps.lonName: if inps.latName and inps.lonName:
geo2radar(geo_file, inps.outfile, inps.latName, inps.lonName) geo2radar(geo_file, inps.outfile, inps.latName, inps.lonName)
return return

15
contrib/stack/stripmapStack/topo.py
View File

@ -401,7 +401,8 @@ def runMultilook(in_dir, out_dir, alks, rlks):
return out_dir return out_dir
def runMultilookGdal(in_dir, out_dir, alks, rlks): def runMultilookGdal(in_dir, out_dir, alks, rlks, in_ext='.rdr', out_ext='.rdr',
fbase_list=['hgt', 'incLocal', 'lat', 'lon', 'los', 'shadowMask', 'waterMask']):
print('generate multilooked geometry files with alks={} and rlks={}'.format(alks, rlks)) print('generate multilooked geometry files with alks={} and rlks={}'.format(alks, rlks))
import gdal import gdal
@ -411,10 +412,9 @@ def runMultilookGdal(in_dir, out_dir, alks, rlks):
print('create directory: {}'.format(out_dir)) print('create directory: {}'.format(out_dir))
# multilook files one by one # multilook files one by one
for fbase in ['hgt', 'incLocal', 'lat', 'lon', 'los', 'shadowMask', 'waterMask']: for fbase in fbase_list:
fname = '{}.rdr'.format(fbase) in_file = os.path.join(in_dir, '{}{}'.format(fbase, in_ext))
in_file = os.path.join(in_dir, fname) out_file = os.path.join(out_dir, '{}{}'.format(fbase, out_ext))
out_file = os.path.join(out_dir, fname)
if os.path.isfile(in_file): if os.path.isfile(in_file):
ds = gdal.Open(in_file, gdal.GA_ReadOnly) ds = gdal.Open(in_file, gdal.GA_ReadOnly)
@ -434,8 +434,9 @@ def runMultilookGdal(in_dir, out_dir, alks, rlks):
# copy the full resolution xml/vrt file from ./merged/geom_master to ./geom_master # copy the full resolution xml/vrt file from ./merged/geom_master to ./geom_master
# to facilitate the number of looks extraction # to facilitate the number of looks extraction
# the file path inside .xml file is not, but should, updated # the file path inside .xml file is not, but should, updated
shutil.copy(in_file+'.xml', out_file+'.full.xml') if in_file != out_file+'.full':
shutil.copy(in_file+'.vrt', out_file+'.full.vrt') shutil.copy(in_file+'.xml', out_file+'.full.xml')
shutil.copy(in_file+'.vrt', out_file+'.full.vrt')
return out_dir return out_dir