Skip to content

Commit

Permalink
Merge pull request #70 from gshiroma/changes_after_calval_point_relea…
Browse files Browse the repository at this point in the history 
…se_6

More changes for the RTC-S1 SAS Release 5 (R5)
  • Loading branch information
@gshiroma
gshiroma authored Sep 9, 2023
2 parents 1b30f97 + 80c2bdd commit 452806a
Show file tree
Hide file tree
Showing 7 changed files with 112 additions and 75 deletions.
2 changes: 1 addition & 1 deletion src/rtc/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -612,7 +612,7 @@ def build_empty_vrt(filename, length, width, fill_value, dtype='Float32',
f' </VRTRasterBand> \n'
f'</VRTDataset> \n')

with open(filename, 'a') as out:
with open(filename, 'w') as out:
out.write(vrt_contents)

if os.path.isfile(filename):
Expand Down
29 changes: 24 additions & 5 deletions src/rtc/h5_prep.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,6 +233,28 @@ def save_hdf5_file(hdf5_obj, output_hdf5_file, clip_max,
logger.info(f'file saved: {output_hdf5_file}')


def get_product_version(product_version_runconfig):
'''
Returns the product version from the product
version defined by the user in the runconfig.
If the runconfig product version is not set, use
the SOFTWARE_VERSION instead
Parameters
---------
product_version_runconfig: scalar
RunConfig product version
Returns
-------
product_version: scalar
Product version
'''
if product_version_runconfig:
return product_version_runconfig
return SOFTWARE_VERSION


def create_hdf5_file(product_id, output_hdf5_file, orbit, burst, cfg,
processing_datetime, is_mosaic):
'''Create HDF5 file
Expand Down Expand Up @@ -375,11 +397,8 @@ def get_metadata_dict(product_id: str,
product_type = cfg_in.groups.primary_executable.product_type

# product version
product_version_float = cfg_in.groups.product_group.product_version
if product_version_float is None:
product_version = SOFTWARE_VERSION
else:
product_version = f'{product_version_float:.1f}'
product_version_runconfig = cfg_in.groups.product_group.product_version
product_version = get_product_version(product_version_runconfig)

# DEM description
dem_file_description = cfg_in.dem_file_description
Expand Down
122 changes: 71 additions & 51 deletions src/rtc/mosaic_geobursts.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ def _compute_distance_to_burst_center(image, geotransform, no_data_value):

x_distance_image, y_distance_image = np.meshgrid(x_vector, y_vector)
distance = np.sqrt((dy * (y_distance_image - center_of_mass[0])) ** 2 +
(dx * (x_distance_image - center_of_mass[1])) ** 2 )
(dx * (x_distance_image - center_of_mass[1])) ** 2)

return distance

Expand All @@ -166,8 +166,9 @@ def is_invalid(input_array, no_data_value):
return np.logical_or(input_array == no_data_value, nan_mask)


def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir='',
geogrid_in=None, temp_files_list=None, verbose=True):
def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode,
scratch_dir='', geogrid_in=None, temp_files_list=None,
verbose=True):
'''
Mosaic S-1 geobursts and return the mosaic as dictionary
Expand All @@ -182,8 +183,9 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
scratch_dir: str (optional)
Directory for temporary files
geogrid_in: isce3.product.GeoGridParameters, default: None
Geogrid information to determine the output mosaic's shape and projection
The geogrid of the output mosaic will automatically determined when it is None
Geogrid information to determine the output mosaic's shape and
projection. If empty, the geogrid of the output mosaic will
determined automatically
temp_files_list: list (optional)
Mutable list of temporary files. If provided,
paths to the temporary files generated will be
Expand Down Expand Up @@ -215,7 +217,7 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
list_geo_transform[i, :] = raster_in.GetGeoTransform()
list_dimension[i, :] = (raster_in.RasterYSize, raster_in.RasterXSize)

# Check if the number of bands are consistent over the input RTC rasters
# Check if the number of bands are consistent with the input rasters
if num_bands is None:
num_bands = raster_in.RasterCount

Expand All @@ -235,15 +237,14 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
no_data_value = None

if geogrid_in is None:
# determine GeoTransformation, posting, dimension, and projection from the input raster
# determine the geogrid from the input raster
for i in range(num_raster):
if list_geo_transform[:, 1].max() == list_geo_transform[:, 1].min():
posting_x = list_geo_transform[0,1]
posting_x = list_geo_transform[0, 1]

if list_geo_transform[:, 5].max() == list_geo_transform[:, 5].min():
posting_y = list_geo_transform[0,5]
posting_y = list_geo_transform[0, 5]

# determine the dimension and the upper left corner of the output mosaic
xmin_mosaic = list_geo_transform[:, 0].min()
ymax_mosaic = list_geo_transform[:, 3].max()
xmax_mosaic = (list_geo_transform[:, 0] +
Expand Down Expand Up @@ -329,7 +330,8 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
gdal_ds = None

# Decide which resampling algorithm to use.
# In case of integer type, use `nearest` e.g. layover shadow mask
# In case of integer type, use `nearest` e.g. layover shadow
# mask
if 'byte' in dtype_name:
resamp_algorithm = 'nearest'
else:
Expand Down Expand Up @@ -357,37 +359,40 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
relocated_file_nlooks = tempfile.NamedTemporaryFile(
dir=scratch_dir, suffix='.tif').name

print(' reprojecting number of looks layer to temporary'
' file:', relocated_file_nlooks)
print(' reprojecting number of looks layer to'
' temporary file:', relocated_file_nlooks)

if temp_files_list is not None:
temp_files_list.append(relocated_file_nlooks)

gdal.Warp(relocated_file_nlooks, path_nlooks,
format='GTiff',
dstSRS=wkt_projection,
outputBounds=[
format='GTiff',
dstSRS=wkt_projection,
outputBounds=[
geogrid_in.start_x,
geogrid_in.start_y +
geogrid_in.length * geogrid_in.spacing_y,
geogrid_in.length * geogrid_in.spacing_y,
geogrid_in.start_x +
geogrid_in.width * geogrid_in.spacing_x,
geogrid_in.width * geogrid_in.spacing_x,
geogrid_in.start_y],
multithread=True,
xRes=geogrid_in.spacing_x,
yRes=abs(geogrid_in.spacing_y),
resampleAlg=resamp_algorithm,
errorThreshold=0,
**warp_options)
multithread=True,
xRes=geogrid_in.spacing_x,
yRes=abs(geogrid_in.spacing_y),
resampleAlg=resamp_algorithm,
errorThreshold=0,
**warp_options)
path_nlooks = relocated_file_nlooks

offset_imgx = 0
offset_imgy = 0
else:

# calculate the burst RTC's offset wrt. the output mosaic in the image coordinate
offset_imgx = int((list_geo_transform[i, 0] - xmin_mosaic) / posting_x + 0.5)
offset_imgy = int((list_geo_transform[i, 3] - ymax_mosaic) / posting_y + 0.5)
# calculate the burst RTC's offset wrt. the output mosaic in the
# image coordinate
offset_imgx = int((list_geo_transform[i, 0] - xmin_mosaic) /
posting_x + 0.5)
offset_imgy = int((list_geo_transform[i, 3] - ymax_mosaic) /
posting_y + 0.5)

# If images start before the mosaic, set up crop start coordinates
if offset_imgx < 0:
Expand All @@ -398,7 +403,8 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
offset_imgy = 0

if verbose:
print(f' image offset (x, y): ({offset_imgx}, {offset_imgy})')
print(f' image offset (x, y): ({offset_imgx},'
f' {offset_imgy})')

if path_nlooks is not None:
nlooks_gdal_ds = gdal.Open(path_nlooks, gdal.GA_ReadOnly)
Expand Down Expand Up @@ -426,28 +432,33 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
if crop_start_x > 0:
arr_rtc = arr_rtc[:, crop_start_x:]

no_data_value = band_ds.GetNoDataValue()

if i == 0 and i_band == 0:

no_data_value = band_ds.GetNoDataValue()
gdal_dtype = band_ds.DataType
dtype_name = gdal.GetDataTypeName(gdal_dtype).lower()
if (IS_MODE_AVERAGE_ENABLED_FOR_BYTE_DTYPE is False and
'byte' in dtype_name and mosaic_mode.lower() == 'average'):
'byte' in dtype_name and mosaic_mode.lower() ==
'average'):
print('WARNING mosaicking of a byte array in the "average"'
' mode is disabled to prevent incorrect averaging of'
' masked values. Using mosaic mode "first" for this'
' layer')
mosaic_mode = 'first'

if mosaic_mode.lower() == 'average':
arr_numerator = np.zeros((num_bands, dim_mosaic[0], dim_mosaic[1]),
arr_numerator = np.zeros((num_bands, dim_mosaic[0],
dim_mosaic[1]),
dtype=arr_rtc.dtype)
arr_denominator = np.zeros(dim_mosaic, dtype=arr_rtc.dtype)
else:
arr_numerator = np.full((num_bands, dim_mosaic[0], dim_mosaic[1]),
arr_numerator = np.full((num_bands, dim_mosaic[0],
dim_mosaic[1]),
no_data_value, dtype=arr_rtc.dtype)
if mosaic_mode.lower() == 'bursts_center':
arr_distance = np.full(dim_mosaic, no_data_value, dtype=arr_rtc.dtype)
arr_distance = np.full(dim_mosaic, no_data_value,
dtype=arr_rtc.dtype)

if i_band == 0:
length = min(arr_rtc.shape[0], dim_mosaic[0] - offset_imgy)
Expand All @@ -464,10 +475,12 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=

if path_nlooks is not None:
arr_denominator[offset_imgy: offset_imgy + length,
offset_imgx: offset_imgx + width] += arr_nlooks
offset_imgx: offset_imgx + width] += \
arr_nlooks
else:
arr_denominator[offset_imgy: offset_imgy + length,
offset_imgx: offset_imgx + width] += np.asarray(
offset_imgx: offset_imgx + width] += \
np.asarray(
np.logical_not(is_invalid(arr_rtc, no_data_value)),
dtype=arr_rtc.dtype)

Expand All @@ -493,18 +506,21 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=
arr_new_distance = _compute_distance_to_burst_center(
arr_rtc, geotransform, no_data_value)

arr_distance_temp = arr_distance[offset_imgy: offset_imgy + length,
offset_imgx: offset_imgx + width]
# update the distance value in `arr_distance` that corresponds to the
# input raster/band, only where corresponding pixel is not valid, and
# the distance array in the input raster is shorter than what is
# in `arr_distance`
ind = np.logical_or(is_invalid(arr_distance_temp, no_data_value),
arr_distance_temp = arr_distance[
offset_imgy: offset_imgy + length,
offset_imgx: offset_imgx + width]
# update the distance value in `arr_distance` that corresponds
# to the input raster/band, only where corresponding pixel is
# not valid, and the distance array in the input raster is
# shorter than what is in `arr_distance`
ind = np.logical_or(is_invalid(arr_distance_temp,
no_data_value),
arr_new_distance <= arr_distance_temp)

arr_distance_temp[ind] = arr_new_distance[ind]
arr_distance[offset_imgy: offset_imgy + length,
offset_imgx: offset_imgx + width] = arr_distance_temp
offset_imgx: offset_imgx + width] = \
arr_distance_temp

del arr_distance_temp

Expand Down Expand Up @@ -544,7 +560,8 @@ def compute_mosaic_array(list_rtc_images, list_nlooks, mosaic_mode, scratch_dir=

def mosaic_single_output_file(list_rtc_images, list_nlooks, mosaic_filename,
mosaic_mode, scratch_dir='', geogrid_in=None,
temp_files_list=None, output_raster_format='GTiff',
temp_files_list=None,
output_raster_format='GTiff',
verbose=True):
'''
Mosaic RTC images saving the output into a single multi-band file
Expand All @@ -560,8 +577,9 @@ def mosaic_single_output_file(list_rtc_images, list_nlooks, mosaic_filename,
scratch_dir: str (optional)
Directory for temporary files
geogrid_in: isce3.product.GeoGridParameters, default: None
Geogrid information to determine the output mosaic's shape and projection
The geogrid of the output mosaic will automatically determined when it is None
Geogrid information to determine the output mosaic's shape and
projection. If empty, the geogrid of the output mosaic will
determined automatically
temp_files_list: list (optional)
Mutable list of temporary files. If provided,
paths to the temporary files generated will be
Expand Down Expand Up @@ -599,7 +617,8 @@ def mosaic_single_output_file(list_rtc_images, list_nlooks, mosaic_filename,
raster_out = drv_out.Create(mosaic_filename, width, length, num_bands,
datatype_mosaic)

raster_out.SetGeoTransform((xmin_mosaic, posting_x, 0, ymax_mosaic, 0, posting_y))
raster_out.SetGeoTransform((xmin_mosaic, posting_x, 0,
ymax_mosaic, 0, posting_y))
raster_out.SetProjection(wkt_projection)

for i_band in range(num_bands):
Expand Down Expand Up @@ -636,8 +655,9 @@ def mosaic_multiple_output_files(
scratch_dir: str (optional)
Directory for temporary files
geogrid_in: isce3.product.GeoGridParameters, default: None
Geogrid information to determine the output mosaic's shape and projection
The geogrid of the output mosaic will automatically determined when it is None
Geogrid information to determine the output mosaic's shape and
projection. If empty, the geogrid of the output mosaic will
determined automatically
temp_files_list: list (optional)
Mutable list of temporary files. If provided,
paths to the temporary files generated will be
Expand Down Expand Up @@ -683,8 +703,8 @@ def mosaic_multiple_output_files(
raster_out = drv_out.Create(output_file, width, length, nbands,
datatype_mosaic)

raster_out.SetGeoTransform((xmin_mosaic, posting_x, 0, ymax_mosaic, 0,
posting_y))
raster_out.SetGeoTransform((xmin_mosaic, posting_x, 0,
ymax_mosaic, 0, posting_y))

raster_out.SetProjection(wkt_projection)

Expand Down
18 changes: 9 additions & 9 deletions src/rtc/rtc_s1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
from rtc.version import VERSION as SOFTWARE_VERSION
from rtc.h5_prep import (save_hdf5_file, create_hdf5_file,
get_metadata_dict,
get_product_version,
all_metadata_dict_to_geotiff_metadata_dict,
layer_hdf5_dict,
DATA_BASE_GROUP,
Expand Down Expand Up @@ -266,13 +267,10 @@ def run_parallel(cfg: RunConfig, logfile_path, flag_logger_full_format):

# primary executable
product_type = cfg.groups.primary_executable.product_type
product_version_float = cfg.groups.product_group.product_version
rtc_s1_static_validity_start_date = \
cfg.groups.product_group.rtc_s1_static_validity_start_date
if product_version_float is None:
product_version = SOFTWARE_VERSION
else:
product_version = f'{product_version_float:.1f}'
product_version_runconfig = cfg.groups.product_group.product_version
product_version = get_product_version(product_version_runconfig)

# unpack processing parameters
processing_namespace = cfg.groups.processing
Expand Down Expand Up @@ -301,7 +299,7 @@ def run_parallel(cfg: RunConfig, logfile_path, flag_logger_full_format):
abs(cfg.geogrid.spacing_y)) / 2)
mosaic_product_id = populate_product_id(
runconfig_product_id, burst_ref, processing_datetime, product_version,
rtc_s1_static_validity_start_date, pixel_spacing_avg, product_type,
pixel_spacing_avg, product_type, rtc_s1_static_validity_start_date,
is_mosaic=True)

# set scratch directory and output_dir
Expand Down Expand Up @@ -375,7 +373,8 @@ def run_parallel(cfg: RunConfig, logfile_path, flag_logger_full_format):
save_incidence_angle = geocode_namespace.save_incidence_angle
save_local_inc_angle = geocode_namespace.save_local_inc_angle
save_projection_angle = geocode_namespace.save_projection_angle
save_rtc_anf_projection_angle = geocode_namespace.save_rtc_anf_projection_angle
save_rtc_anf_projection_angle = \
geocode_namespace.save_rtc_anf_projection_angle
save_range_slope = geocode_namespace.save_range_slope
save_nlooks = geocode_namespace.save_nlooks

Expand Down Expand Up @@ -811,7 +810,7 @@ def run_parallel(cfg: RunConfig, logfile_path, flag_logger_full_format):
if save_rtc_anf_gamma0_to_sigma0:
output_metadata_dict[
LAYER_NAME_RTC_ANF_GAMMA0_TO_SIGMA0][1].append(
rtc_anf_file)
rtc_anf_gamma0_to_sigma0_file)

radar_grid_file_dict = {}

Expand All @@ -821,7 +820,8 @@ def run_parallel(cfg: RunConfig, logfile_path, flag_logger_full_format):
LAYER_NAME_INCIDENCE_ANGLE: save_incidence_angle,
LAYER_NAME_LOCAL_INCIDENCE_ANGLE: save_local_inc_angle,
LAYER_NAME_PROJECTION_ANGLE: save_projection_angle,
LAYER_NAME_RTC_ANF_PROJECTION_ANGLE: save_rtc_anf_projection_angle,
LAYER_NAME_RTC_ANF_PROJECTION_ANGLE:
save_rtc_anf_projection_angle,
LAYER_NAME_RANGE_SLOPE: save_range_slope,
LAYER_NAME_DEM: save_dem}

Expand Down
Loading

0 comments on commit 452806a

Please sign in to comment.