Skip to content

Commit

Permalink
Move make_geotiff to library
Browse files Browse the repository at this point in the history
  • Loading branch information
@yumorishita
yumorishita committed Feb 9, 2021
1 parent b102c74 commit b39223f
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 85 deletions.
46 changes: 36 additions & 10 deletions LiCSBAS_lib/LiCSBAS_io_lib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -24,6 +26,7 @@
import subprocess as subp
import datetime as dt
import statsmodels.api as sm
from osgeo import gdal, osr


#%%
Expand All @@ -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:
Expand Down Expand Up @@ -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 = ''
Expand Down Expand Up @@ -106,15 +132,15 @@ 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
Return: bperp
"""
bperp = []
bperp_dict = {}

### Determine type of bperp_file; old or not
with open(bperp_file) as f:
line = f.readline().split() #list
Expand All @@ -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


Expand All @@ -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


Expand All @@ -171,7 +197,7 @@ def read_ifg_list(ifg_listfile):
ifgdates.append(ifgd)
line = f.readline()
f.close()

return ifgdates


Expand Down
82 changes: 31 additions & 51 deletions bin/LiCSBAS_decomposeLOS.py
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -43,55 +47,31 @@
import numpy as np
import time
from decimal import Decimal
import LiCSBAS_io_lib as io_lib

class Usage(Exception):
"""Usage context manager"""
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)


#%% 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']
Expand Down Expand Up @@ -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... ')
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down
40 changes: 16 additions & 24 deletions bin/LiCSBAS_flt2geotiff.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python3
"""
v1.5 20200902 Yu Morishita, GSI
v1.5.1 20210209 Yu Morishita, GSI
========
Overview
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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')
Expand All @@ -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
Expand Down

0 comments on commit b39223f

Please sign in to comment.