From b39223f431cbb29fd53c3f8d3aef1952a4898055 Mon Sep 17 00:00:00 2001 From: Yu Morishita Date: Tue, 9 Feb 2021 10:12:57 +0900 Subject: [PATCH] Move make_geotiff to library --- LiCSBAS_lib/LiCSBAS_io_lib.py | 46 +++++++++++++++----- bin/LiCSBAS_decomposeLOS.py | 82 +++++++++++++---------------------- bin/LiCSBAS_flt2geotiff.py | 40 +++++++---------- 3 files changed, 83 insertions(+), 85 deletions(-) diff --git a/LiCSBAS_lib/LiCSBAS_io_lib.py b/LiCSBAS_lib/LiCSBAS_io_lib.py index 3e78844..e133e0c 100644 --- a/LiCSBAS_lib/LiCSBAS_io_lib.py +++ b/LiCSBAS_lib/LiCSBAS_io_lib.py @@ -8,6 +8,8 @@ ========= Changelog ========= +v1.3 20210209 Yu Morioshita, GSI + - Add make_geotiff v1.2.1 20201211 Yu Morioshita, GSI - Skip invalid lines in baselines file in read_bperp_file v1.2 20200703 Yu Morioshita, GSI @@ -24,6 +26,7 @@ import subprocess as subp import datetime as dt import statsmodels.api as sm +from osgeo import gdal, osr #%% @@ -37,10 +40,33 @@ def make_dummy_bperp(bperp_file, imdates): elif np.mod(i, 4)==0: bp = -np.random.rand()/2 #-0.5~0 ifg_dt = dt.datetime.strptime(imd, '%Y%m%d').toordinal() - dt.datetime.strptime(imdates[0], '%Y%m%d').toordinal() - + print('{:3d} {} {} {:5.2f} {:4d} {} {:4d} {} {:5.2f}'.format(i, imdates[0], imd, bp, ifg_dt, 0, ifg_dt, 0, bp), file=f) +#%% +def make_geotiff(data, latn_p, lonw_p, dlat, dlon, outfile, compress_option): + length, width = data.shape + if data.dtype == np.float32: + dtype = gdal.GDT_Float32 + nodata = np.nan ## or 0? + elif data.dtype == np.uint8: + dtype = gdal.GDT_Byte + nodata = None + + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(outfile, width, length, 1, dtype, options=compress_option) + outRaster.SetGeoTransform((lonw_p, dlon, 0, latn_p, 0, dlat)) + outband = outRaster.GetRasterBand(1) + outband.WriteArray(data) + if nodata is not None: outband.SetNoDataValue(nodata) + outRaster.SetMetadataItem('AREA_OR_POINT', 'Area') + outRasterSRS = osr.SpatialReference() + outRasterSRS.ImportFromEPSG(4326) + outRaster.SetProjection(outRasterSRS.ExportToWkt()) + outband.FlushCache() + + #%% def make_point_kml(lat, lon, kmlfile): with open(kmlfile, "w") as f: @@ -69,7 +95,7 @@ def make_tstxt(x, y, imdates, ts, tsfile, refx1, refx2, refy1, refy2, gap, lat=N imdates_yr = (imdates_ordinal-imdates_ordinal[0])/365.25 A = sm.add_constant(imdates_yr) #[1, t] vconst, vel = sm.OLS(ts, A, missing='drop').fit().params - + ### Identify gaps ixs_gap = np.where(gap==1)[0] # n_im-1, bool gap_str = '' @@ -106,7 +132,7 @@ def read_bperp_file(bperp_file, imdates): Old bperp_file contains (m: primary (master), s:secondary, sm: single prime): - num mdate sdate bp dt dt_m_sm dt_s_sm bp_m_sm bp_s_sm + num mdate sdate bp dt dt_m_sm dt_s_sm bp_m_sm bp_s_sm 1 20170218 20170326 96.6 36.0 -12.0 24.0 34.2 130.9 2 20170302 20170314 32.4 12.0 0.0 12.0 0.0 32.4 @@ -114,7 +140,7 @@ def read_bperp_file(bperp_file, imdates): """ bperp = [] bperp_dict = {} - + ### Determine type of bperp_file; old or not with open(bperp_file) as f: line = f.readline().split() #list @@ -125,20 +151,20 @@ def read_bperp_file(bperp_file, imdates): for l in f: if len(l.split()) == 4: bperp_dict[l.split()[1]] = l.split()[2] - + else: ## old format with open(bperp_file) as f: for l in f: bperp_dict[l.split()[1]] = l.split()[-2] bperp_dict[l.split()[2]] = l.split()[-1] - + for imd in imdates: if imd in bperp_dict: bperp.append(float(bperp_dict[imd])) else: ## If no key exists print('ERROR: bperp for {} not found!'.format(imd), file=sys.stderr) return False - + return bperp @@ -148,12 +174,12 @@ def read_img(file, length, width, dtype=np.float32, endian='little'): Read image data into numpy array. endian: 'little' or 'big' (not 'little' is regarded as 'big') """ - + if endian == 'little': data = np.fromfile(file, dtype=dtype).reshape((length, width)) else: data = np.fromfile(file, dtype=dtype).byteswap().reshape((length, width)) - + return data @@ -171,7 +197,7 @@ def read_ifg_list(ifg_listfile): ifgdates.append(ifgd) line = f.readline() f.close() - + return ifgdates diff --git a/bin/LiCSBAS_decomposeLOS.py b/bin/LiCSBAS_decomposeLOS.py index ba99664..ca0517b 100755 --- a/bin/LiCSBAS_decomposeLOS.py +++ b/bin/LiCSBAS_decomposeLOS.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -v1.1.1 20201112 Yu Morishita, GSI +v1.1.2 20210209 Yu Morishita, GSI This script decomposes 2 (or more) LOS displacement data to EW and UD components assuming no NS displacement (neglecting NS). Positive values in the decomposed data mean eastward and upward displacement. The multiple LOS input data can have different coverage and resolution as they are resampled to the common area and resolution during the processing. @@ -13,12 +13,16 @@ ===== LiCSBAS_decomposeLOS.py -f files.txt [-o outfile] [-r resampleAlg] [--out_stats] - -f Text file containing input GeoTIFF file paths of LOS displacement - (or velocity), E and N components of LOS unit vector + -f Text file containing input GeoTIFF file paths of LOS displacement + (or velocity) and E and N components of LOS unit vector from >=2 directions + Note: GeoTIFF files can be created by LiCSBAS_flt2geotiff.py Format: dispfile1 Efile1 Nfile1 dispfile2 Efile2 Nfile2 ... + Example: + 046D_vel.geo.tif 046D_E.geo.tif 046D_N.geo.tif + 068A_vel.geo.tif 068A_E.geo.tif 068A_N.geo.tif -o Prefix of output decomposed file (Default: no prefix, [EW|UD].geo.tif) -r Resampling algorithm (Default: bilinear) (see https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r) @@ -27,8 +31,8 @@ """ #%% Change log ''' -v1.1.1 20201112 Yu Morishita, GSI - - Small bug fix +v1.1.2 20210209 Yu Morishita, GSI + - Move make_geotiff to library v1.1 20200608 Yu Morishita, GSI - Add --out_stats option v1.0 20200528 Yu Morishita, GSI @@ -43,6 +47,7 @@ import numpy as np import time from decimal import Decimal +import LiCSBAS_io_lib as io_lib class Usage(Exception): """Usage context manager""" @@ -50,40 +55,15 @@ def __init__(self, msg): self.msg = msg -#%% -def make_geotiff(data, length, width, latn_p, lonw_p, dlat, dlon, outfile, compress_option): - - if data.dtype == np.float32: - dtype = gdal.GDT_Float32 - nodata = np.nan ## or 0? - elif data.dtype == np.uint8: - dtype = gdal.GDT_Byte - nodata = None - - driver = gdal.GetDriverByName('GTiff') - outRaster = driver.Create(outfile, width, length, 1, dtype, options=compress_option) - outRaster.SetGeoTransform((lonw_p, dlon, 0, latn_p, 0, dlat)) - outband = outRaster.GetRasterBand(1) - outband.WriteArray(data) - if nodata is not None: outband.SetNoDataValue(nodata) - outRaster.SetMetadataItem('AREA_OR_POINT', 'Point') - outRasterSRS = osr.SpatialReference() - outRasterSRS.ImportFromEPSG(4326) - outRaster.SetProjection(outRasterSRS.ExportToWkt()) - outband.FlushCache() - - return - - #%% Main def main(argv=None): - + #%% Check argv if argv == None: argv = sys.argv - + start = time.time() - ver='1.1.1'; date=20201112; author="Y. Morishita" + ver='1.1.2'; date=20210209; author="Y. Morishita" print("\n{} ver{} {} {}".format(os.path.basename(argv[0]), ver, date, author), flush=True) print("{} {}".format(os.path.basename(argv[0]), ' '.join(argv[1:])), flush=True) @@ -91,7 +71,7 @@ def main(argv=None): #%% Set default infiletxt = [] out_prefix = '' - resampleAlg = 'bilinear' #'cubicspline'# 'near' # 'cubic' + resampleAlg = 'bilinear' #'cubicspline'# 'near' # 'cubic' out_stats_flag = False compress_option = ['COMPRESS=DEFLATE', 'PREDICTOR=3'] compress_option_uint = ['COMPRESS=DEFLATE', 'PREDICTOR=1'] @@ -144,7 +124,7 @@ def main(argv=None): n_data = len(data_tifs) print('\nNumber of input LOS data: {}'.format(n_data)) - + #%% Identify area with at least 1 each from E and W ### All latlon values in this script are in pixel registration print('\nRead area of each GeoTIFF... ') @@ -203,7 +183,7 @@ def main(argv=None): lon_e = lon_e_E if lon_e_E < lon_e_W else lon_e_W lat_s = lat_s_E if lat_s_E > lat_s_W else lat_s_W lat_n = lat_n_E if lat_n_E < lat_n_W else lat_n_W - + width = int((lon_e-lon_w)/dlon) lon_e = lon_w + dlon*width length = int((lat_s-lat_n)/dlat) @@ -276,10 +256,10 @@ def main(argv=None): det = (a11*a22-a12**2) det[det==0] = np.nan ## To avoid zero division detinv = 1/det - + ew_part = detinv*(a22*be-a12*bu) ud_part = detinv*(-a12*be+a11*bu) - + ew = np.zeros_like(bool_valid, dtype=np.float32)*np.nan ew[bool_valid] = ew_part ud = np.zeros_like(bool_valid, dtype=np.float32)*np.nan @@ -289,31 +269,31 @@ def main(argv=None): #%% Save geotiff outfileEW = out_prefix + 'EW.geo.tif' outfileUD = out_prefix + 'UD.geo.tif' - make_geotiff(ew, length, width, lat_n, lon_w, dlat, dlon, outfileEW, compress_option) - make_geotiff(ud, length, width, lat_n, lon_w, dlat, dlon, outfileUD, compress_option) + io_lib.make_geotiff(ew, lat_n, lon_w, dlat, dlon, outfileEW, compress_option) + io_lib.make_geotiff(ud, lat_n, lon_w, dlat, dlon, outfileUD, compress_option) #%% Stats if out_stats_flag: if n_data >= 3: for i in range(n_data): - outfile_resid = out_prefix + 'resid_LOS{}.geo.tif'.format(i+1) + outfile_resid = out_prefix + 'resid_LOS{}.geo.tif'.format(i+1) data_part_list[i][data_part_list[i]==0] = np.nan resid_los_part = data_part_list[i] - \ (LOSe_part_list[i]*ew_part + LOSu_part_list[i]*ud_part) resid_los = np.zeros_like(bool_valid, dtype=np.float32)*np.nan - resid_los[bool_valid] = resid_los_part - make_geotiff(resid_los, length, width, lat_n, lon_w, dlat, dlon, outfile_resid, compress_option) + resid_los[bool_valid] = resid_los_part + io_lib.make_geotiff(resid_los, lat_n, lon_w, dlat, dlon, outfile_resid, compress_option) ### n_data - outfile_n_data = out_prefix + 'n_data_fromE.geo.tif' - make_geotiff(n_data_fromE, length, width, lat_n, lon_w, dlat, dlon, outfile_n_data, compress_option_uint) - outfile_n_data = out_prefix + 'n_data_fromW.geo.tif' - make_geotiff(n_data_fromW, length, width, lat_n, lon_w, dlat, dlon, outfile_n_data, compress_option_uint) - outfile_n_data = out_prefix + 'n_data_total.geo.tif' - make_geotiff(n_data_total, length, width, lat_n, lon_w, dlat, dlon, outfile_n_data, compress_option_uint) - - + outfile_n_data = out_prefix + 'n_data_fromE.geo.tif' + io_lib.make_geotiff(n_data_fromE, lat_n, lon_w, dlat, dlon, outfile_n_data, compress_option_uint) + outfile_n_data = out_prefix + 'n_data_fromW.geo.tif' + io_lib.make_geotiff(n_data_fromW, lat_n, lon_w, dlat, dlon, outfile_n_data, compress_option_uint) + outfile_n_data = out_prefix + 'n_data_total.geo.tif' + io_lib.make_geotiff(n_data_total, lat_n, lon_w, dlat, dlon, outfile_n_data, compress_option_uint) + + #%% Finish elapsed_time = time.time()-start hour = int(elapsed_time/3600) diff --git a/bin/LiCSBAS_flt2geotiff.py b/bin/LiCSBAS_flt2geotiff.py index 3aff0a8..f9e3f24 100755 --- a/bin/LiCSBAS_flt2geotiff.py +++ b/bin/LiCSBAS_flt2geotiff.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 """ -v1.5 20200902 Yu Morishita, GSI +v1.5.1 20210209 Yu Morishita, GSI ======== Overview @@ -26,6 +26,8 @@ #%% Change log ''' +v1.5.1 20210209 Yu Morishita, GSI + - Move make_geotiff to library v1.5 20200902 Yu Morishita, GSI - Do not add .geo when already added v1.4 20200214 Yu Morishita, Uni of Leeds and GSI @@ -58,13 +60,13 @@ def __init__(self, msg): #%% Main def main(argv=None): - + #%% Check argv if argv == None: argv = sys.argv - + start = time.time() - ver=1.5; date=20200902; author="Y. Morishita" + ver="1.5.1"; date=20210209; author="Y. Morishita" print("\n{} ver{} {} {}".format(os.path.basename(argv[0]), ver, date, author), flush=True) print("{} {}".format(os.path.basename(argv[0]), ' '.join(argv[1:])), flush=True) @@ -153,16 +155,16 @@ def main(argv=None): data.byteswap().tofile(infile+'.bigendian') infile = infile+'.bigendian' - + print('{} {} {} {} {}'.format(gamma_comm, dempar, infile, dtype, outfile)) call=[gamma_comm, dempar, infile, dtype, outfile] subp.run(call) - + if endian == 'little': ### Remove temporary file os.remove(infile) - - + + #%% if gdal and osr available else: print('Use gdal module') @@ -174,35 +176,25 @@ def main(argv=None): width = int(io_lib.get_param_par(dempar, 'width')) length = int(io_lib.get_param_par(dempar, 'nlines')) - + dlat = float(io_lib.get_param_par(dempar, 'post_lat')) dlon = float(io_lib.get_param_par(dempar, 'post_lon')) - + lat_n_g = float(io_lib.get_param_par(dempar, 'corner_lat')) #grid reg lon_w_g = float(io_lib.get_param_par(dempar, 'corner_lon')) #grid reg - + ## Grid registration to pixel registration by shifing half pixel lat_n_p = lat_n_g - dlat/2 lon_w_p = lon_w_g - dlon/2 - + data = io_lib.read_img(infile, length, width, endian=endian) if zero2nan_flag: ### Replace 0 with nan data[data==0] = np.nan if nan2zero_flag: ### Replace nan with 0 data[np.isnan(data)] = 0 - - driver = gdal.GetDriverByName('GTiff') - outRaster = driver.Create(outfile, width, length, 1, gdal.GDT_Float32, options=compress_option) - outRaster.SetGeoTransform((lon_w_p, dlon, 0, lat_n_p, 0, dlat)) - outband = outRaster.GetRasterBand(1) - outband.WriteArray(data) - if nodata is not None: outband.SetNoDataValue(nodata) - outRaster.SetMetadataItem('AREA_OR_POINT', 'Point') - outRasterSRS = osr.SpatialReference() - outRasterSRS.ImportFromEPSG(4326) - outRaster.SetProjection(outRasterSRS.ExportToWkt()) - outband.FlushCache() + + io_lib.make_geotiff(data, lat_n_p, lon_w_p, dlat, dlon, outfile, compress_option) #%% Finish