From 639e33e6142912e6385c372ba1bbd779abefc862 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Mon, 7 Aug 2023 16:18:56 -0700 Subject: [PATCH 01/11] Test routine and main framework of the code --- src/compass/mosaic_cslc_backscatter.py | 253 +++++++++++++++++++++++++ 1 file changed, 253 insertions(+) create mode 100644 src/compass/mosaic_cslc_backscatter.py diff --git a/src/compass/mosaic_cslc_backscatter.py b/src/compass/mosaic_cslc_backscatter.py new file mode 100644 index 00000000..b7720079 --- /dev/null +++ b/src/compass/mosaic_cslc_backscatter.py @@ -0,0 +1,253 @@ +''' +Script to mosaic the input CSLC rasters +''' +import os + +import numpy as np +from osgeo import gdal, osr +from collections import Counter +from itertools import repeat + +PATH_CSLC_LAYER_IN_HDF = '/data' +PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' + +def get_most_frequent_epsg(gdal_raster_list: list): + ''' + Get the most frequent EPSG among the rasters in 'gdal_raster_list' + + Parameters + ---------- + gdal_raster_list: list + List of strings for GDAL raster dataset + + Returns + ------- + most_common_epsg: int + Most common EPSG code among the input parameters + + ''' + # Get a list of all EPSG codes + epsg_list = [_get_epsg(raster) for raster in gdal_raster_list] + + # Count the occurrence of each EPSG code + counter = Counter(epsg_list) + + # Return the most common EPSG code + most_common_epsg = counter.most_common(1)[0][0] + + return most_common_epsg + + +def _get_epsg(gdal_raster_path: str): + ''' + Get the EPSG of the input raster + ''' + ds_input = gdal.Open(gdal_raster_path, gdal.GA_ReadOnly) + projection = ds_input.GetProjection() + srs = osr.SpatialReference(wkt=projection) + return int(srs.GetAuthorityCode(None)) + + +def get_cslc_gdal_dataset(cslc_path, pol, epsg_out, dx, dy, snap=0, + noise_correction=True, + resampling_alg='bicubic'): + ''' + Get the GDAL dataset for CSLC layer and resampled & reprojected noise LUT (when user opted in) + + ''' + # Get the geotransform and dimensions + prefix_netcdf = f'DERIVED_SUBDATASET:AMPLITUDE:NETCDF' + + path_cslc = f'{prefix_netcdf}:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' + path_noise_lut = f'NETCDF:{cslc_path}:{PATH_NOISE_LAYER_IN_HDF}/thermal_noise_lut' + + ds_in = gdal.Open(path_cslc, gdal.GA_ReadOnly) + epsg_cslc = _get_epsg(path_cslc) + + ds_noise_lut = (gdal.Open(path_noise_lut, gdal.GA_ReadOnly) if + noise_correction else None) + + if epsg_cslc == epsg_out: + # Do not perform reprojection, + # and return the Dataset object diretly from the original + ds_cslc_layer = ds_in + return (ds_cslc_layer, ds_noise_lut) + + # Reproject the CSLC layer and + # Define the source and target spatial references + gt_in = ds_in.GetGeoTransform() + xsize = ds_in.RasterXSize + ysize = ds_in.RasterYSize + + src_srs = osr.SpatialReference() + src_srs.ImportFromWkt(ds_in.GetProjectionRef()) + dst_srs = osr.SpatialReference() + dst_srs.ImportFromEPSG(epsg_out) + + # Define a coordinate transformation + transform = osr.CoordinateTransformation(src_srs, dst_srs) + + # Define the corner points + corners = [(0, 0), (0, ysize), (xsize, ysize), (xsize, 0)] + + # Convert the corner points in `epsg_out` + transformed_corners = [] + for corner in corners: + x, y = (gt_in[0] + corner[0] * gt_in[1], + gt_in[3] + corner[1] * gt_in[5]) + x_transformed, y_transformed, _ = transform.TransformPoint(x, y) + transformed_corners.append([x_transformed, y_transformed]) + transformed_corners = np.array(transformed_corners) + + if snap != 0: + transformed_corners = np.round(transformed_corners / snap) * snap + + # Extract the extent of the reprojected raster + extent_reprojected = [np.min(transformed_corners[:, 0]), + np.min(transformed_corners[:, 1]), + np.max(transformed_corners[:, 0]), + np.max(transformed_corners[:, 1])] + + # Put together the warp options + warp_options = gdal.WarpOptions(format='MEM', + resampleAlg=resampling_alg, + xRes=dx, + yRes=abs(dy), + outputBounds=extent_reprojected, + dstSRS=dst_srs) + + ds_cslc_reprojected = gdal.Warp('', ds_in, options=warp_options) + ds_noise_reprojected = (gdal.Warp('', ds_noise_lut, options=warp_options) if + ds_noise_lut else None) + + # Return the warped dataset + return (ds_cslc_reprojected, ds_noise_reprojected) + + +def apply_noise_correction(): + ''' + Apply noise correction using the geocoded noist LUT in the product + ''' + pass + + +def mosaic_list_order_mode(gdal_raster_list, mosaic_path): + pass + + +def mosaic_nearest_centroid_mode(gdal_raster_list, mosaic_path): + pass + + +def compute_mosaic_geotransform_dimension(ds_list): + """ + Compute the GeoTransform vector that covers all rasters in `raster_list` + """ + + # Initialize list to store extents + if len(ds_list) == 0: + raise RuntimeError('Empty ds_list was provided') + + extents = [] + for ds in ds_list: + gt = ds.GetGeoTransform() + # Compute extent from GeoTransform + xmin = gt[0] + xmax = gt[0] + gt[1] * ds.RasterXSize + ymin = gt[3] + gt[5] * ds.RasterYSize + ymax = gt[3] + extents.append([xmin, xmax, ymin, ymax]) + + extents = np.array(extents) + + # Find common extent + xmin = np.min(extents[:, 0]) + xmax = np.max(extents[:, 1]) + ymin = np.min(extents[:, 2]) + ymax = np.max(extents[:, 3]) + + # Create and return GeoTransform + common_gt = (xmin, gt[1], gt[2], ymax, gt[4], gt[5]) + + width_mosaic = int((xmax - xmin)/ gt[1] + 0.5) + height_mosaic = int((ymin - ymax)/ gt[5] + 0.5) + return common_gt, (width_mosaic, height_mosaic) + + + +def run(cslc_path_list, pol, mosaic_path, + dx_mosaic=None, dy_mosaic=None, snap_meters=30, + epsg_mosaic=None, mode='list_order', + apply_noise_correcion=True, resampling_alg='bicubic'): + ''' + Workflow of the CSLC backscatter mosaic + + Parameters + ---------- + cslc_path_list + pol: str + Polarization to mosaic + mosaic_path: + Path to the output mosaic + dx_mosaic: float, Optional + Spacing of the mosaic in x direction + dy_mosaic: float, Optional + Spacing of the mosaic in y direction + snap_meters: + Snapping in meters + epsg_mosaic: int, optional + EPSG of the output mosaic. + If None, then the EPSG will be determined based on the input bursts + mode: + Mosaic mode + apply_noise_correcion: bool + Flag whether or not to apply thermal noise correction + resampling_alg: str + Resampling algorithm, if necessary + + ''' + + # Determine the EPSG of the mosaic when it is not provided + if epsg_mosaic is None: + epsg_mosaic = get_most_frequent_epsg( + [f'NETCDF:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' for cslc_path in cslc_path_list]) + + if dx_mosaic is None or dy_mosaic is None: + ds_1st_cslc = gdal.Open(f'NETCDF:{cslc_path_list[0]}:{PATH_CSLC_LAYER_IN_HDF}/{pol}') + gt_1st_cslc = ds_1st_cslc.GetGeoTransform() + dx_mosaic = gt_1st_cslc[1] if dx_mosaic is None else dx_mosaic + dy_mosaic = gt_1st_cslc[5] if dy_mosaic is None else dy_mosaic + + # Get the list of gdal raster dataset + num_raster_in = len(cslc_path_list) + cslc_amp_ds_list = [None] * num_raster_in + noise_lut_ds_list = [None] * num_raster_in + + for i_cslc, cslc_path in enumerate(cslc_path_list): + print(f'Loading CSLC layer in amplutude, and nouise LUTs: {i_cslc + 1} / {len(cslc_path_list)}', end='\r') + datasets_cslc = get_cslc_gdal_dataset(cslc_path, pol, epsg_mosaic, + dx_mosaic, dy_mosaic, snap_meters, + apply_noise_correcion, resampling_alg) + cslc_amp_ds_list[i_cslc] = datasets_cslc[0] + noise_lut_ds_list[i_cslc] = datasets_cslc[1] + + print(' ') + # Get the size of mosaic + gt_mosaic, (width_mosaic, height_mosaic) = compute_mosaic_geotransform_dimension(cslc_amp_ds_list) + + print('sdfsadf') + + # Compute the array for the mosaicked image + + + + # Write out the mosaicked image + + +if __name__=='__main__': + import glob + HOMEDIR = os.getenv('HOME') + burst_cslc_dir = os.path.join(HOMEDIR, 'Documents/OPERA_SCRATCH/CSLC/MOSAIC_TEST_SITE/output_s1_cslc','**/**/t*.h5') + output_mosaic_path = os.path.join(HOMEDIR, 'Documents/OPERA_SCRATCH/CSLC/MOSAIC_TEST_SITE/output_s1_cslc/mosaic.tif') + list_burst = glob.glob(burst_cslc_dir) + run(list_burst, 'VV', output_mosaic_path) From fd9fd8a8ebcad387fd3566de4d69f79e9d09ba39 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 8 Aug 2023 18:00:42 -0700 Subject: [PATCH 02/11] ongoing CLI implementation; applying radiometric normalization; --- src/compass/mosaic_cslc_backscatter.py | 347 ++++++++++++++++++++----- 1 file changed, 280 insertions(+), 67 deletions(-) diff --git a/src/compass/mosaic_cslc_backscatter.py b/src/compass/mosaic_cslc_backscatter.py index b7720079..c8746614 100644 --- a/src/compass/mosaic_cslc_backscatter.py +++ b/src/compass/mosaic_cslc_backscatter.py @@ -1,16 +1,59 @@ ''' Script to mosaic the input CSLC rasters ''' +import argparse import os +from collections import Counter import numpy as np from osgeo import gdal, osr -from collections import Counter -from itertools import repeat +# Constants for dataset location in HDF5 file PATH_CSLC_LAYER_IN_HDF = '/data' +PATH_LOCAL_INCIDENCE_ANGLE = '/data/local_incidence_angle' PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' + +def get_parser(): + ''' + Get the parser for CLI + ''' + parser = argparse.ArgumentParser( + description='Comparison script with burst ID', + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('-c', + dest='cslc_list_file', + type=str, + default='', + help='cslc list file') + + parser.add_argument('-s', + dest='static_list_file', + type=str, + default=[], + help='static list file') + + parser.add_argument('-o', + dest='mosaic_out_path', + type=str, + default=None, + help='path to the burst DB') + + parser.add_argument('-cd', + dest='cslc_directory', + type=str, + default=None, + help='path to the cslc directory') + + parser.add_argument('-sd', + dest='cslc_static_directory', + type=str, + default=None, + help='path to the cslc static layer directory') + + return parser + + def get_most_frequent_epsg(gdal_raster_list: list): ''' Get the most frequent EPSG among the rasters in 'gdal_raster_list' @@ -48,17 +91,19 @@ def _get_epsg(gdal_raster_path: str): return int(srs.GetAuthorityCode(None)) -def get_cslc_gdal_dataset(cslc_path, pol, epsg_out, dx, dy, snap=0, +def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, snap=0, noise_correction=True, - resampling_alg='bicubic'): + radiometric_normalization=True, + resampling_alg='BILINEAR'): ''' Get the GDAL dataset for CSLC layer and resampled & reprojected noise LUT (when user opted in) ''' # Get the geotransform and dimensions - prefix_netcdf = f'DERIVED_SUBDATASET:AMPLITUDE:NETCDF' + prefix_netcdf = 'DERIVED_SUBDATASET:AMPLITUDE:NETCDF' path_cslc = f'{prefix_netcdf}:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' + path_local_incidence = f'NETCDF:{cslc_static_path}:{PATH_LOCAL_INCIDENCE_ANGLE}' path_noise_lut = f'NETCDF:{cslc_path}:{PATH_NOISE_LAYER_IN_HDF}/thermal_noise_lut' ds_in = gdal.Open(path_cslc, gdal.GA_ReadOnly) @@ -67,13 +112,11 @@ def get_cslc_gdal_dataset(cslc_path, pol, epsg_out, dx, dy, snap=0, ds_noise_lut = (gdal.Open(path_noise_lut, gdal.GA_ReadOnly) if noise_correction else None) - if epsg_cslc == epsg_out: - # Do not perform reprojection, - # and return the Dataset object diretly from the original - ds_cslc_layer = ds_in - return (ds_cslc_layer, ds_noise_lut) + ds_local_incidence_angle = (gdal.Open(path_local_incidence, gdal.GA_ReadOnly) if + radiometric_normalization else None) - # Reproject the CSLC layer and + # Prepare for the reprojection + # Reproject the CSLC layer, and # Define the source and target spatial references gt_in = ds_in.GetGeoTransform() xsize = ds_in.RasterXSize @@ -92,21 +135,29 @@ def get_cslc_gdal_dataset(cslc_path, pol, epsg_out, dx, dy, snap=0, # Convert the corner points in `epsg_out` transformed_corners = [] - for corner in corners: - x, y = (gt_in[0] + corner[0] * gt_in[1], - gt_in[3] + corner[1] * gt_in[5]) - x_transformed, y_transformed, _ = transform.TransformPoint(x, y) - transformed_corners.append([x_transformed, y_transformed]) - transformed_corners = np.array(transformed_corners) - - if snap != 0: - transformed_corners = np.round(transformed_corners / snap) * snap - - # Extract the extent of the reprojected raster - extent_reprojected = [np.min(transformed_corners[:, 0]), - np.min(transformed_corners[:, 1]), - np.max(transformed_corners[:, 0]), - np.max(transformed_corners[:, 1])] + + # Compute the extet of the resampled raster + if epsg_cslc == epsg_out: + extent_reprojected=[gt_in[0], + gt_in[3] + gt_in[5] * ysize, + gt_in[0] + gt_in[1] * xsize, + gt_in[3]] + else: + for corner in corners: + x_from, y_from = (gt_in[0] + corner[0] * gt_in[1], + gt_in[3] + corner[1] * gt_in[5]) + x_transformed, y_transformed, _ = transform.TransformPoint(x_from, y_from) + transformed_corners.append([x_transformed, y_transformed]) + transformed_corners = np.array(transformed_corners) + + if snap != 0: + transformed_corners = np.round(transformed_corners / snap) * snap + + # Extract the extent of the reprojected raster + extent_reprojected = [np.min(transformed_corners[:, 0]), + np.min(transformed_corners[:, 1]), + np.max(transformed_corners[:, 0]), + np.max(transformed_corners[:, 1])] # Put together the warp options warp_options = gdal.WarpOptions(format='MEM', @@ -115,47 +166,146 @@ def get_cslc_gdal_dataset(cslc_path, pol, epsg_out, dx, dy, snap=0, yRes=abs(dy), outputBounds=extent_reprojected, dstSRS=dst_srs) + # Resample the noise LUT + ds_noise_resampled = (gdal.Warp('', ds_noise_lut, options=warp_options) if + ds_noise_lut else None) + + if epsg_cslc == epsg_out: + # Do not perform reprojection when EPSG of the burst is the same as that of the mosaic + # return the GDAL Dataset objects as they are, except for noise LUT + return (ds_in, ds_local_incidence_angle, ds_noise_resampled) + ds_incidence_angle_reprojected = (gdal.Warp('', ds_local_incidence_angle, options=warp_options) if ds_local_incidence_angle else None) ds_cslc_reprojected = gdal.Warp('', ds_in, options=warp_options) - ds_noise_reprojected = (gdal.Warp('', ds_noise_lut, options=warp_options) if - ds_noise_lut else None) + # Return the warped dataset - return (ds_cslc_reprojected, ds_noise_reprojected) + return (ds_cslc_reprojected, ds_incidence_angle_reprojected, ds_noise_resampled) -def apply_noise_correction(): +def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): ''' Apply noise correction using the geocoded noist LUT in the product ''' - pass + # Load CSLC amplutude into array + arr_cslc = ds_cslc_amp.ReadAsArray() + + if ds_noise_lut is not None: + print('Applying noise removal') + # Resample noise LUT to the geogrid of the CSLC amplitude + #geotransform_cslc = ds_cslc_amp.GetGeoTransform() + #width_cslc = ds_cslc_amp.RasterXSize + #height_cslc = ds_cslc_amp.RasterYSize + + #bbox_cslc=[geotransform_cslc[0], + # geotransform_cslc[3] + geotransform_cslc[5] * height_cslc, + # geotransform_cslc[0] + geotransform_cslc[1] * width_cslc, + # geotransform_cslc[3]] + + #ds_resampled_lut = gdal.Warp( + # destNameOrDestDS='', + # srcDSOrSrcDSTab=ds_noise_lut, + # format='MEM', + # outputBounds=bbox_cslc, + # width=width_cslc, + # height=height_cslc) + + # Apply noise correction + + arr_cslc = arr_cslc ** 2 - ds_noise_lut.ReadAsArray() + arr_cslc[arr_cslc<0.0] = 0.0 + arr_cslc = np.sqrt(arr_cslc) + if ds_local_incidence_angle is not None: + print('Applying radiometric mormalization') + # Apply radiometric normalization using local incidence angle + correction_factor = np.cos(np.deg2rad(ds_local_incidence_angle.ReadAsArray())) + arr_cslc /= correction_factor -def mosaic_list_order_mode(gdal_raster_list, mosaic_path): - pass + return arr_cslc + + +def find_upperleft_pixel_index(geotransform_burst, geotransform_mosaic): + ''' + Find the upperleft corner of the burst image in mosaic image, in pixels + + Parameters + ---------- + geotransform_burst: tuple + Geotransform parameter for burst + geotransform_mosaic: tuple + Geotransform parameter for mosaic + + Returns + ------- + upperleft_y_px, upperleft_x_px: int + relative location of upper-left corner of the burst in mosaic + ''' + xmin_burst = geotransform_burst[0] + xmin_mosaic = geotransform_mosaic[0] + + ymax_burst = geotransform_burst[3] + ymax_mosaic = geotransform_mosaic[3] + + upperleft_x_px = int((xmin_burst - xmin_mosaic) / geotransform_mosaic[1]) + upperleft_y_px = int((ymax_burst - ymax_mosaic) / geotransform_mosaic[5]) + + return (upperleft_y_px, upperleft_x_px) + + +def mosaic_list_order_mode(cslc_raster_list, local_incidence_angle_list, noise_lut_raster_list, + geotransform_mosaic, shape_mosaic): + ''' + Perform the mosaicking in "list order mode" + i.e. The pixels from the earlier raster gets replaced by the following rasters + ''' + array_mosaic = np.zeros(shape_mosaic, dtype=np.float32) / 0 -def mosaic_nearest_centroid_mode(gdal_raster_list, mosaic_path): - pass + for i_raster, cslc_raster in enumerate(cslc_raster_list): + print(f'Processing: {i_raster + 1} of {len(cslc_raster_list)}') + noise_lut_raster = noise_lut_raster_list[i_raster] + local_incidence_angle_raster = local_incidence_angle_list[i_raster] + uly, ulx = find_upperleft_pixel_index(cslc_raster.GetGeoTransform(), geotransform_mosaic) + amplitude_burst = load_amplitude(cslc_raster, local_incidence_angle_raster, noise_lut_raster) + length_burst, width_burst = amplitude_burst.shape + subset_array_mosaic = array_mosaic[uly:uly+length_burst, ulx:ulx+width_burst] + mask_valid = ~np.isnan(amplitude_burst) + subset_array_mosaic[mask_valid] = amplitude_burst[mask_valid] -def compute_mosaic_geotransform_dimension(ds_list): + return array_mosaic + + +def mosaic_nearest_centroid_mode(cslc_raster_list, noise_lut_raster_list, + geotransform_mosaic, shape_mosaic): + ''' + Perform the mosaicking in "nearest centroid" mode + + ''' + array_mosaic = np.array(shape_mosaic, dtype=np.float32) + + + return array_mosaic + + +def compute_mosaic_geotransform_dimension(ds_burst_list): """ Compute the GeoTransform vector that covers all rasters in `raster_list` """ # Initialize list to store extents - if len(ds_list) == 0: - raise RuntimeError('Empty ds_list was provided') + if len(ds_burst_list) == 0: + raise RuntimeError('Empty ds_burst_list was provided') extents = [] - for ds in ds_list: - gt = ds.GetGeoTransform() + for ds_burst in ds_burst_list: + gt_burst = ds_burst.GetGeoTransform() # Compute extent from GeoTransform - xmin = gt[0] - xmax = gt[0] + gt[1] * ds.RasterXSize - ymin = gt[3] + gt[5] * ds.RasterYSize - ymax = gt[3] + xmin = gt_burst[0] + xmax = gt_burst[0] + gt_burst[1] * ds_burst.RasterXSize + ymin = gt_burst[3] + gt_burst[5] * ds_burst.RasterYSize + ymax = gt_burst[3] extents.append([xmin, xmax, ymin, ymax]) extents = np.array(extents) @@ -167,18 +317,54 @@ def compute_mosaic_geotransform_dimension(ds_list): ymax = np.max(extents[:, 3]) # Create and return GeoTransform - common_gt = (xmin, gt[1], gt[2], ymax, gt[4], gt[5]) + gt_burst = ds_burst_list[0].GetGeoTransform() + gt_mosaic = (xmin, gt_burst[1], gt_burst[2], ymax, gt_burst[4], gt_burst[5]) + + width_mosaic = int((xmax - xmin)/ gt_burst[1] + 0.5) + height_mosaic = int((ymin - ymax)/ gt_burst[5] + 0.5) + return gt_mosaic, (width_mosaic, height_mosaic) + + +def save_mosaicked_array(mosaic_arr, geotransform_mosaic, epsg_mosaic, mosaic_filename): + ''' + Save mosaicked array into GDAL Raster + + Parameters + ---------- + mosaic_arr: np.ndarray + Mosaicked image as numpy array + geotransform_mosaic: tuple + Geotransform parameter for mosaic + epsg_mosaic: + EPSG for mosaic as projection parameter + mosaic_filename: + File name of the mosaic raster to save + ''' - width_mosaic = int((xmax - xmin)/ gt[1] + 0.5) - height_mosaic = int((ymin - ymax)/ gt[5] + 0.5) - return common_gt, (width_mosaic, height_mosaic) + # Create the projection from the EPSG code + srs_mosaic = osr.SpatialReference() + srs_mosaic.ImportFromEPSG(epsg_mosaic) + projection_mosaic = srs_mosaic.ExportToWkt() + length_mosaic, width_mosaic = mosaic_arr.shape + driver = gdal.GetDriverByName('GTiff') + ds_out = driver.Create(mosaic_filename, + width_mosaic, length_mosaic, 1, gdal.GDT_Float32, + options=['COMPRESS=LZW', 'BIGTIFF=YES']) + ds_out.SetGeoTransform(geotransform_mosaic) + ds_out.SetProjection(projection_mosaic) -def run(cslc_path_list, pol, mosaic_path, + ds_out.GetRasterBand(1).WriteArray(mosaic_arr) + + ds_out = None + + +def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, dx_mosaic=None, dy_mosaic=None, snap_meters=30, - epsg_mosaic=None, mode='list_order', - apply_noise_correcion=True, resampling_alg='bicubic'): + epsg_mosaic=None, mosaic_mode='list_order', + apply_noise_correcion=True, apply_radiometric_normalization=True, + resampling_alg='bicubic'): ''' Workflow of the CSLC backscatter mosaic @@ -198,19 +384,20 @@ def run(cslc_path_list, pol, mosaic_path, epsg_mosaic: int, optional EPSG of the output mosaic. If None, then the EPSG will be determined based on the input bursts - mode: + mosaic_mode: Mosaic mode apply_noise_correcion: bool Flag whether or not to apply thermal noise correction resampling_alg: str Resampling algorithm, if necessary - ''' # Determine the EPSG of the mosaic when it is not provided if epsg_mosaic is None: + print('EPSG was not provided. Finding most common projection among the input rasters.') epsg_mosaic = get_most_frequent_epsg( [f'NETCDF:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' for cslc_path in cslc_path_list]) + print(f'EPSG of the mosaic will be: {epsg_mosaic}') if dx_mosaic is None or dy_mosaic is None: ds_1st_cslc = gdal.Open(f'NETCDF:{cslc_path_list[0]}:{PATH_CSLC_LAYER_IN_HDF}/{pol}') @@ -222,32 +409,58 @@ def run(cslc_path_list, pol, mosaic_path, num_raster_in = len(cslc_path_list) cslc_amp_ds_list = [None] * num_raster_in noise_lut_ds_list = [None] * num_raster_in - + local_incidence_angle_ds_list = [None] * num_raster_in + for i_cslc, cslc_path in enumerate(cslc_path_list): - print(f'Loading CSLC layer in amplutude, and nouise LUTs: {i_cslc + 1} / {len(cslc_path_list)}', end='\r') - datasets_cslc = get_cslc_gdal_dataset(cslc_path, pol, epsg_mosaic, + cslc_static_path = cslc_static_path_list[i_cslc] + print('Loading CSLC layer in amplutude, and nouise LUTs: ' + f'{i_cslc + 1} / {len(cslc_path_list)}', + end='\r') + datasets_cslc = get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_mosaic, dx_mosaic, dy_mosaic, snap_meters, - apply_noise_correcion, resampling_alg) + apply_noise_correcion, apply_radiometric_normalization, + resampling_alg) cslc_amp_ds_list[i_cslc] = datasets_cslc[0] - noise_lut_ds_list[i_cslc] = datasets_cslc[1] + local_incidence_angle_ds_list[i_cslc] = datasets_cslc[1] + noise_lut_ds_list[i_cslc] = datasets_cslc[2] print(' ') # Get the size of mosaic - gt_mosaic, (width_mosaic, height_mosaic) = compute_mosaic_geotransform_dimension(cslc_amp_ds_list) - - print('sdfsadf') + gt_mosaic, (width_mosaic, height_mosaic) = \ + compute_mosaic_geotransform_dimension(cslc_amp_ds_list) # Compute the array for the mosaicked image + mosaic_functions = { + 'list_order': mosaic_list_order_mode, + 'nearest': mosaic_nearest_centroid_mode + } + if mosaic_mode not in mosaic_functions: + raise NotImplementedError(f'Mosaicking was not implemented for mode: {mosaic_mode}') + mosaic_arr = mosaic_functions[mosaic_mode]( + cslc_amp_ds_list, local_incidence_angle_ds_list, noise_lut_ds_list, gt_mosaic, + (height_mosaic, width_mosaic)) # Write out the mosaicked image + print(f'Writing mosaic to: {mosaic_path}') + save_mosaicked_array(mosaic_arr, gt_mosaic, epsg_mosaic, mosaic_path) + + +def main(): + ''' + Entrypoint of the fiunction + ''' + parser = get_parser() + args = parser.parse_args() + + run(list_burst, list_burst_static, 'VV', output_mosaic_path) + + + + if __name__=='__main__': - import glob - HOMEDIR = os.getenv('HOME') - burst_cslc_dir = os.path.join(HOMEDIR, 'Documents/OPERA_SCRATCH/CSLC/MOSAIC_TEST_SITE/output_s1_cslc','**/**/t*.h5') - output_mosaic_path = os.path.join(HOMEDIR, 'Documents/OPERA_SCRATCH/CSLC/MOSAIC_TEST_SITE/output_s1_cslc/mosaic.tif') - list_burst = glob.glob(burst_cslc_dir) - run(list_burst, 'VV', output_mosaic_path) + main() + From c7aa54d2c37ef9437899de3dc5392b51cb2b0481 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 8 Aug 2023 21:48:16 -0700 Subject: [PATCH 03/11] relocation of the mosaicking script; CLI implementation --- .../{ => utils}/mosaic_cslc_backscatter.py | 187 ++++++++++-------- 1 file changed, 108 insertions(+), 79 deletions(-) rename src/compass/{ => utils}/mosaic_cslc_backscatter.py (80%) mode change 100644 => 100755 diff --git a/src/compass/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py old mode 100644 new mode 100755 similarity index 80% rename from src/compass/mosaic_cslc_backscatter.py rename to src/compass/utils/mosaic_cslc_backscatter.py index c8746614..0a1e08fe --- a/src/compass/mosaic_cslc_backscatter.py +++ b/src/compass/utils/mosaic_cslc_backscatter.py @@ -1,8 +1,9 @@ +#! python + ''' Script to mosaic the input CSLC rasters ''' import argparse -import os from collections import Counter import numpy as np @@ -14,46 +15,6 @@ PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' -def get_parser(): - ''' - Get the parser for CLI - ''' - parser = argparse.ArgumentParser( - description='Comparison script with burst ID', - formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('-c', - dest='cslc_list_file', - type=str, - default='', - help='cslc list file') - - parser.add_argument('-s', - dest='static_list_file', - type=str, - default=[], - help='static list file') - - parser.add_argument('-o', - dest='mosaic_out_path', - type=str, - default=None, - help='path to the burst DB') - - parser.add_argument('-cd', - dest='cslc_directory', - type=str, - default=None, - help='path to the cslc directory') - - parser.add_argument('-sd', - dest='cslc_static_directory', - type=str, - default=None, - help='path to the cslc static layer directory') - - return parser - - def get_most_frequent_epsg(gdal_raster_list: list): ''' Get the most frequent EPSG among the rasters in 'gdal_raster_list' @@ -175,12 +136,15 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, sn # return the GDAL Dataset objects as they are, except for noise LUT return (ds_in, ds_local_incidence_angle, ds_noise_resampled) - ds_incidence_angle_reprojected = (gdal.Warp('', ds_local_incidence_angle, options=warp_options) if ds_local_incidence_angle else None) + ds_incidence_angle_reprojected = \ + (gdal.Warp('', ds_local_incidence_angle, options=warp_options) + if ds_local_incidence_angle else None) ds_cslc_reprojected = gdal.Warp('', ds_in, options=warp_options) - # Return the warped dataset - return (ds_cslc_reprojected, ds_incidence_angle_reprojected, ds_noise_resampled) + return (ds_cslc_reprojected, + ds_incidence_angle_reprojected, + ds_noise_resampled) def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): @@ -192,24 +156,6 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None if ds_noise_lut is not None: print('Applying noise removal') - # Resample noise LUT to the geogrid of the CSLC amplitude - #geotransform_cslc = ds_cslc_amp.GetGeoTransform() - #width_cslc = ds_cslc_amp.RasterXSize - #height_cslc = ds_cslc_amp.RasterYSize - - #bbox_cslc=[geotransform_cslc[0], - # geotransform_cslc[3] + geotransform_cslc[5] * height_cslc, - # geotransform_cslc[0] + geotransform_cslc[1] * width_cslc, - # geotransform_cslc[3]] - - #ds_resampled_lut = gdal.Warp( - # destNameOrDestDS='', - # srcDSOrSrcDSTab=ds_noise_lut, - # format='MEM', - # outputBounds=bbox_cslc, - # width=width_cslc, - # height=height_cslc) - # Apply noise correction arr_cslc = arr_cslc ** 2 - ds_noise_lut.ReadAsArray() @@ -219,7 +165,8 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None if ds_local_incidence_angle is not None: print('Applying radiometric mormalization') # Apply radiometric normalization using local incidence angle - correction_factor = np.cos(np.deg2rad(ds_local_incidence_angle.ReadAsArray())) + correction_factor = \ + np.cos(np.deg2rad(ds_local_incidence_angle.ReadAsArray())) arr_cslc /= correction_factor return arr_cslc @@ -253,8 +200,11 @@ def find_upperleft_pixel_index(geotransform_burst, geotransform_mosaic): return (upperleft_y_px, upperleft_x_px) -def mosaic_list_order_mode(cslc_raster_list, local_incidence_angle_list, noise_lut_raster_list, - geotransform_mosaic, shape_mosaic): +def mosaic_list_order_mode(cslc_raster_list, + local_incidence_angle_list, + noise_lut_raster_list, + geotransform_mosaic, + shape_mosaic): ''' Perform the mosaicking in "list order mode" i.e. The pixels from the earlier raster gets replaced by the following rasters @@ -267,7 +217,9 @@ def mosaic_list_order_mode(cslc_raster_list, local_incidence_angle_list, noise_l noise_lut_raster = noise_lut_raster_list[i_raster] local_incidence_angle_raster = local_incidence_angle_list[i_raster] uly, ulx = find_upperleft_pixel_index(cslc_raster.GetGeoTransform(), geotransform_mosaic) - amplitude_burst = load_amplitude(cslc_raster, local_incidence_angle_raster, noise_lut_raster) + amplitude_burst = load_amplitude(cslc_raster, + local_incidence_angle_raster, + noise_lut_raster) length_burst, width_burst = amplitude_burst.shape subset_array_mosaic = array_mosaic[uly:uly+length_burst, ulx:ulx+width_burst] @@ -277,16 +229,14 @@ def mosaic_list_order_mode(cslc_raster_list, local_incidence_angle_list, noise_l return array_mosaic -def mosaic_nearest_centroid_mode(cslc_raster_list, noise_lut_raster_list, +def mosaic_nearest_centroid_mode(cslc_raster_list, local_incidence_angle_list, + noise_lut_raster_list, geotransform_mosaic, shape_mosaic): ''' Perform the mosaicking in "nearest centroid" mode ''' - array_mosaic = np.array(shape_mosaic, dtype=np.float32) - - - return array_mosaic + raise NotImplementedError('This mosaic function was not implemented yet') def compute_mosaic_geotransform_dimension(ds_burst_list): @@ -364,7 +314,7 @@ def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, dx_mosaic=None, dy_mosaic=None, snap_meters=30, epsg_mosaic=None, mosaic_mode='list_order', apply_noise_correcion=True, apply_radiometric_normalization=True, - resampling_alg='bicubic'): + resampling_alg='BILINEAR'): ''' Workflow of the CSLC backscatter mosaic @@ -416,13 +366,15 @@ def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, print('Loading CSLC layer in amplutude, and nouise LUTs: ' f'{i_cslc + 1} / {len(cslc_path_list)}', end='\r') - datasets_cslc = get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_mosaic, + datasets_cslc = get_cslc_gdal_dataset(cslc_path, cslc_static_path, + pol, epsg_mosaic, dx_mosaic, dy_mosaic, snap_meters, - apply_noise_correcion, apply_radiometric_normalization, + apply_noise_correcion, + apply_radiometric_normalization, resampling_alg) - cslc_amp_ds_list[i_cslc] = datasets_cslc[0] - local_incidence_angle_ds_list[i_cslc] = datasets_cslc[1] - noise_lut_ds_list[i_cslc] = datasets_cslc[2] + (cslc_amp_ds_list[i_cslc], + local_incidence_angle_ds_list[i_cslc], + noise_lut_ds_list[i_cslc]) = datasets_cslc print(' ') # Get the size of mosaic @@ -447,6 +399,75 @@ def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, save_mosaicked_array(mosaic_arr, gt_mosaic, epsg_mosaic, mosaic_path) + +def get_parser(): + ''' + Get the parser for CLI + ''' + parser = argparse.ArgumentParser( + description='Comparison script with burst ID', + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('-c', + dest='cslc_list_file', + type=str, + default='', + help='cslc list file') + + parser.add_argument('-s', + dest='cslc_static_list_file', + type=str, + default=[], + help='static list file') + + parser.add_argument('-o', + dest='mosaic_out_path', + type=str, + default=None, + help='path to the burst DB') + + parser.add_argument('-p', + dest='pol', + type=str, + default='VV', + help='Polarization') + + parser.add_argument('-te', + dest='target_spacing', + type=float, + nargs=2, + default=[None, None], + help=('spacing of the mosaic. If not provided, ' + 'the spacing of the first input CSLC will be used.'), + metavar=('NUM1', 'NUM2')) + + parser.add_argument('--snap', + dest='snap', + type=float, + default=30.0, + help='Snapping value of the mosaic grid.') + + parser.add_argument('--mode', + dest='mosaic_mode', + type=str, + default='list_order', + help='Mosaic mode') + + parser.add_argument('--noise_correction_off', + dest='apply_noise_correction', + default=True, + action='store_false', + help='Turn off the noise correction') + + parser.add_argument('--radiometric_normalization_off', + dest='apply_radiometric_normalization', + default=True, + action='store_false', + help='Turn off the radiometric normalization') + + + return parser + + def main(): ''' Entrypoint of the fiunction @@ -454,7 +475,16 @@ def main(): parser = get_parser() args = parser.parse_args() - run(list_burst, list_burst_static, 'VV', output_mosaic_path) + with (open(args.cslc_list_file, 'r') as fin_cslc, + open(args.cslc_static_list_file, 'r') as fin_cslc_static): + list_cslc = fin_cslc.read().rstrip('\n').split('\n') + list_cslc_static = fin_cslc_static.read().rstrip('\n').split('\n') + + run(list_cslc, list_cslc_static, args.pol, args.mosaic_out_path, + dx_mosaic=args.target_spacing[0], dy_mosaic=args.target_spacing[1], + snap_meters=args.snap, mosaic_mode=args.mosaic_mode, + apply_noise_correcion=args.apply_noise_correction, + apply_radiometric_normalization=args.apply_radiometric_normalization) @@ -463,4 +493,3 @@ def main(): if __name__=='__main__': main() - From 6ab235ecc0cdeb1164c7c7d96f425de7f0ea0455 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Tue, 8 Aug 2023 22:16:55 -0700 Subject: [PATCH 04/11] Docstrings added --- src/compass/utils/mosaic_cslc_backscatter.py | 116 +++++++++++++++---- 1 file changed, 96 insertions(+), 20 deletions(-) diff --git a/src/compass/utils/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py index 0a1e08fe..2d9be323 100755 --- a/src/compass/utils/mosaic_cslc_backscatter.py +++ b/src/compass/utils/mosaic_cslc_backscatter.py @@ -53,14 +53,47 @@ def _get_epsg(gdal_raster_path: str): def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, snap=0, - noise_correction=True, - radiometric_normalization=True, - resampling_alg='BILINEAR'): + noise_correction=True, + radiometric_normalization=True, + resampling_alg='BILINEAR'): ''' - Get the GDAL dataset for CSLC layer and resampled & reprojected noise LUT (when user opted in) + Get the GDAL dataset for CSLC layer and resample & + reprojecte noise LUT (when user opted in) + Parameters + ---------- + cslc_path: str + path to the CSLC HDF5 file + cslc_static_path: str + path to the CSLC static layer HDF5 file + pol: str + Polarization + epsg_out: int + EPSG of the outputs + dx, dy: float + x / y spacings of the outputs + snap: float + Snapping value when the input rasters needs to be resampled + noise_correction: bool + Flag to turn on/off noise correction + radiometric_normalization:bool + Flag to turn on/ott radiometric normalization + resampling_alg: str + Resampling algorithm for gdal.Warp() + + Returns + ------- + (ds_cslc, + ds_incidence_angle, + ds_noise_resampled): tuple + Respectively, GDAL raster dataset for CSLC, + GDAL raster dataset for local incidence angle, and + GDAL raster dataset for noise LUT + (resampled to the same geogrid as the other two) ''' # Get the geotransform and dimensions + + prefix_netcdf = 'DERIVED_SUBDATASET:AMPLITUDE:NETCDF' path_cslc = f'{prefix_netcdf}:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' @@ -88,6 +121,11 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, sn dst_srs = osr.SpatialReference() dst_srs.ImportFromEPSG(epsg_out) + # decide whether or not to warp CSLC and static layer + flag_warp_cslc_and_static = ((epsg_cslc != epsg_out) + or (gt_in[0] != dx) + or (abs(gt_in[3]) != abs(dy))) + # Define a coordinate transformation transform = osr.CoordinateTransformation(src_srs, dst_srs) @@ -131,8 +169,10 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, sn ds_noise_resampled = (gdal.Warp('', ds_noise_lut, options=warp_options) if ds_noise_lut else None) - if epsg_cslc == epsg_out: - # Do not perform reprojection when EPSG of the burst is the same as that of the mosaic + if not flag_warp_cslc_and_static: + # Do not perform reprojection when + # EPSG of the burst is the same as that of the mosaic AND + # the spacing of the CSLC is the same as the mosaic's. # return the GDAL Dataset objects as they are, except for noise LUT return (ds_in, ds_local_incidence_angle, ds_noise_resampled) @@ -149,8 +189,25 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, sn def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): ''' - Apply noise correction using the geocoded noist LUT in the product + Load the amplitude from CSLC GDAL dataset. + Apply noise correction and/or radiometric normalization when + the data are provided. + + Parameters + ---------- + ds_cslc_amp: osgeo.gdal.Dataset + GDAL Raster dataset for CSLC amplitude + ds_local_incidence_angle: osgeo.gdal.Dataset + GDAL Raster dataset for local incidence angle + ds_noise_lut: osgeo.gdal.Dataset + GDAL Raster dataset noise LUT + + Returns + ------- + arr_cslc: np.ndarray + CSLC amplitude ''' + # Load CSLC amplutude into array arr_cslc = ds_cslc_amp.ReadAsArray() @@ -208,6 +265,24 @@ def mosaic_list_order_mode(cslc_raster_list, ''' Perform the mosaicking in "list order mode" i.e. The pixels from the earlier raster gets replaced by the following rasters + + Parameters + ---------- + cslc_raster_list: list + List of GDAL raster dataset for CSLC + local_incidence_angle_list: list + List of GDAL raster dataset for local incidence angle + noise_lut_raster_list: list + List of GDAL raster dataset for noise LUT + geotransform_mosaic: tuple + Geotramsform parameter for the output mosaic + shape_mosaic: tuple + Shape of the output mosaic array + + Returns + ------- + array_mosaic: np.ndarray + Output mosaic as numpy array ''' array_mosaic = np.zeros(shape_mosaic, dtype=np.float32) / 0 @@ -216,13 +291,16 @@ def mosaic_list_order_mode(cslc_raster_list, print(f'Processing: {i_raster + 1} of {len(cslc_raster_list)}') noise_lut_raster = noise_lut_raster_list[i_raster] local_incidence_angle_raster = local_incidence_angle_list[i_raster] - uly, ulx = find_upperleft_pixel_index(cslc_raster.GetGeoTransform(), geotransform_mosaic) + uly, ulx = find_upperleft_pixel_index(cslc_raster.GetGeoTransform(), + geotransform_mosaic) amplitude_burst = load_amplitude(cslc_raster, local_incidence_angle_raster, noise_lut_raster) length_burst, width_burst = amplitude_burst.shape - subset_array_mosaic = array_mosaic[uly:uly+length_burst, ulx:ulx+width_burst] + subset_array_mosaic = array_mosaic[uly : uly + length_burst, + ulx : ulx + width_burst] + mask_valid = ~np.isnan(amplitude_burst) subset_array_mosaic[mask_valid] = amplitude_burst[mask_valid] @@ -240,10 +318,9 @@ def mosaic_nearest_centroid_mode(cslc_raster_list, local_incidence_angle_list, def compute_mosaic_geotransform_dimension(ds_burst_list): - """ - Compute the GeoTransform vector that covers all rasters in `raster_list` - """ - + '''Compute the GeoTransform vector that covers all rasters in `ds_burst_list` + ''' + # Initialize list to store extents if len(ds_burst_list) == 0: raise RuntimeError('Empty ds_burst_list was provided') @@ -275,7 +352,8 @@ def compute_mosaic_geotransform_dimension(ds_burst_list): return gt_mosaic, (width_mosaic, height_mosaic) -def save_mosaicked_array(mosaic_arr, geotransform_mosaic, epsg_mosaic, mosaic_filename): +def save_mosaicked_array(mosaic_arr, geotransform_mosaic, epsg_mosaic, + mosaic_filename): ''' Save mosaicked array into GDAL Raster @@ -344,9 +422,11 @@ def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, # Determine the EPSG of the mosaic when it is not provided if epsg_mosaic is None: - print('EPSG was not provided. Finding most common projection among the input rasters.') + print('EPSG was not provided. ' + 'Finding most common projection among the input rasters.') epsg_mosaic = get_most_frequent_epsg( - [f'NETCDF:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' for cslc_path in cslc_path_list]) + [f'NETCDF:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' + for cslc_path in cslc_path_list]) print(f'EPSG of the mosaic will be: {epsg_mosaic}') if dx_mosaic is None or dy_mosaic is None: @@ -487,9 +567,5 @@ def main(): apply_radiometric_normalization=args.apply_radiometric_normalization) - - - - if __name__=='__main__': main() From 17242bb09af70cc012f41b870e513b61fe53bc02 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Thu, 10 Aug 2023 11:48:42 -0700 Subject: [PATCH 05/11] script for amplitude extraction / correction only; changing the equation for radiometric norrmalization --- src/compass/utils/correct_cslc_amplitude.py | 260 +++++++++++++++++++ src/compass/utils/mosaic_cslc_backscatter.py | 9 +- 2 files changed, 265 insertions(+), 4 deletions(-) create mode 100755 src/compass/utils/correct_cslc_amplitude.py diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py new file mode 100755 index 00000000..5e683b8c --- /dev/null +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -0,0 +1,260 @@ +#! python + +''' +Script to mosaic the input CSLC rasters +''' +import argparse + +import numpy as np +from osgeo import gdal, osr + +# Constants for dataset location in HDF5 file +PATH_CSLC_LAYER_IN_HDF = '/data' +PATH_LOCAL_INCIDENCE_ANGLE = '/data/local_incidence_angle' +PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' + + +def _get_epsg(gdal_raster_path: str): + ''' + Get the EPSG of the input raster + ''' + ds_input = gdal.Open(gdal_raster_path, gdal.GA_ReadOnly) + projection = ds_input.GetProjection() + srs = osr.SpatialReference(wkt=projection) + return int(srs.GetAuthorityCode(None)) + + +def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): + ''' + Load the amplitude from CSLC GDAL dataset. + Apply noise correction and/or radiometric normalization when + the data are provided. + + Parameters + ---------- + ds_cslc_amp: osgeo.gdal.Dataset + GDAL Raster dataset for CSLC amplitude + ds_local_incidence_angle: osgeo.gdal.Dataset + GDAL Raster dataset for local incidence angle + ds_noise_lut: osgeo.gdal.Dataset + GDAL Raster dataset noise LUT + + Returns + ------- + arr_cslc: np.ndarray + CSLC amplitude + ''' + + # Load CSLC amplutude into array + arr_cslc = ds_cslc_amp.ReadAsArray() + + if ds_noise_lut is not None: + print('Applying noise removal') + # Apply noise correction + + arr_cslc = arr_cslc ** 2 - ds_noise_lut.ReadAsArray() + arr_cslc[arr_cslc<0.0] = 0.0 + arr_cslc = np.sqrt(arr_cslc) + + if ds_local_incidence_angle is not None: + print('Applying radiometric mormalization') + local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) + # Apply radiometric normalization using cotangent(local incidence angle) + correction_factor = (np.sin(local_incidence_angle_arr_rad) + / np.cos(local_incidence_angle_arr_rad)) + arr_cslc *= correction_factor + + return arr_cslc + + +def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, + noise_correction=True, + radiometric_normalization=True, + resampling_alg='BILINEAR'): + ''' + Get the GDAL dataset for CSLC layer and resample & + reprojecte noise LUT (when user opted in) + + Parameters + ---------- + cslc_path: str + path to the CSLC HDF5 file + cslc_static_path: str + path to the CSLC static layer HDF5 file + pol: str + Polarization + noise_correction: bool + Flag to turn on/off noise correction + radiometric_normalization:bool + Flag to turn on/ott radiometric normalization + resampling_alg: str + Resampling algorithm for gdal.Warp() + + Returns + ------- + (ds_cslc, + ds_incidence_angle, + ds_noise_resampled): tuple + Respectively, GDAL raster dataset for CSLC, + GDAL raster dataset for local incidence angle, and + GDAL raster dataset for noise LUT + (resampled to the same geogrid as the other two) + ''' + # Get the geotransform and dimensions + + + prefix_netcdf = 'DERIVED_SUBDATASET:AMPLITUDE:NETCDF' + + path_cslc = f'{prefix_netcdf}:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' + path_local_incidence = f'NETCDF:{cslc_static_path}:{PATH_LOCAL_INCIDENCE_ANGLE}' + path_noise_lut = f'NETCDF:{cslc_path}:{PATH_NOISE_LAYER_IN_HDF}/thermal_noise_lut' + + ds_in = gdal.Open(path_cslc, gdal.GA_ReadOnly) + + ds_noise_lut = (gdal.Open(path_noise_lut, gdal.GA_ReadOnly) if + noise_correction else None) + + ds_local_incidence_angle = (gdal.Open(path_local_incidence, gdal.GA_ReadOnly) if + radiometric_normalization else None) + + # Prepare for the reprojection + # Reproject the CSLC layer, and + # Define the source and target spatial references + gt_in = ds_in.GetGeoTransform() + xsize = ds_in.RasterXSize + ysize = ds_in.RasterYSize + + src_srs = osr.SpatialReference() + src_srs.ImportFromWkt(ds_in.GetProjectionRef()) + + # Convert the corner points in `epsg_out` + extent_reprojected=[gt_in[0], + gt_in[3] + gt_in[5] * ysize, + gt_in[0] + gt_in[1] * xsize, + gt_in[3]] + + # Put together the warp options + warp_options = gdal.WarpOptions(format='MEM', + resampleAlg=resampling_alg, + xRes=gt_in[1], + yRes=abs(gt_in[5]), + outputBounds=extent_reprojected, + dstSRS=src_srs) + + # Resample the noise LUT + ds_noise_resampled = (gdal.Warp('', ds_noise_lut, options=warp_options) if + ds_noise_lut else None) + + return (ds_in, + ds_local_incidence_angle, + ds_noise_resampled) + + +def save_amplitude(amplitude_arr, geotransform_cslc, epsg_cslc, + amplitude_filename): + ''' + Save mosaicked array into GDAL Raster + + Parameters + ---------- + amplitude_arr: np.ndarray + Mosaicked image as numpy array + geotransform_cslc: tuple + Geotransform parameter for mosaic + epsg_cslc: + EPSG for mosaic as projection parameter + amplitude_filename: + File name of the mosaic raster to save + ''' + + # Create the projection from the EPSG code + srs_out = osr.SpatialReference() + srs_out.ImportFromEPSG(epsg_cslc) + projection_out = srs_out.ExportToWkt() + + length_mosaic, width_mosaic = amplitude_arr.shape + driver = gdal.GetDriverByName('GTiff') + ds_out = driver.Create(amplitude_filename, + width_mosaic, length_mosaic, 1, gdal.GDT_Float32, + options=['COMPRESS=LZW', 'BIGTIFF=YES']) + + ds_out.SetGeoTransform(geotransform_cslc) + ds_out.SetProjection(projection_out) + + ds_out.GetRasterBand(1).WriteArray(amplitude_arr) + + ds_out = None + + +def get_parser(): + ''' + Get the parser for CLI + ''' + parser = argparse.ArgumentParser( + description='Comparison script with burst ID', + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('-c', + dest='cslc_path', + type=str, + default='', + help='cslc file') + + parser.add_argument('-s', + dest='cslc_static_path', + type=str, + default=[], + help='CSLC static layer file') + + parser.add_argument('-o', + dest='out_path', + type=str, + default=None, + help='path to the output file') + + parser.add_argument('-p', + dest='pol', + type=str, + default='VV', + help='Polarization') + + parser.add_argument('--noise_correction_off', + dest='apply_noise_correction', + default=True, + action='store_false', + help='Turn off the noise correction') + + parser.add_argument('--radiometric_normalization_off', + dest='apply_radiometric_normalization', + default=True, + action='store_false', + help='Turn off the radiometric normalization') + + return parser + + +def main(): + ''' + Entrypoint of the fiunction + ''' + parser = get_parser() + args = parser.parse_args() + + cslc_layer, noise_lut, local_incidence_angle = \ + get_cslc_gdal_dataset(args.cslc_path, + args.cslc_static_path, + args.pol, + args.apply_noise_correction, + args.apply_radiometric_normalization) + amplitude_arr = load_amplitude(cslc_layer, noise_lut, local_incidence_angle) + + geotransform_cslc = cslc_layer.GetGeoTransform() + + proj_cslc = cslc_layer.GetProjection() + srs_cslc = osr.SpatialReference(wkt=proj_cslc) + epsg_cslc = int(srs_cslc.GetAuthorityCode(None)) + + save_amplitude(amplitude_arr, geotransform_cslc, epsg_cslc, args.out_path) + + +if __name__=='__main__': + main() diff --git a/src/compass/utils/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py index 2d9be323..40235701 100755 --- a/src/compass/utils/mosaic_cslc_backscatter.py +++ b/src/compass/utils/mosaic_cslc_backscatter.py @@ -221,10 +221,11 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None if ds_local_incidence_angle is not None: print('Applying radiometric mormalization') - # Apply radiometric normalization using local incidence angle - correction_factor = \ - np.cos(np.deg2rad(ds_local_incidence_angle.ReadAsArray())) - arr_cslc /= correction_factor + local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) + # Apply radiometric normalization using cotangent(local incidence angle) + correction_factor = (np.sin(local_incidence_angle_arr_rad) + / np.cos(local_incidence_angle_arr_rad)) + arr_cslc *= correction_factor return arr_cslc From 8a63510dac92f418f2bbaa152aa4fdb963fce339 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Thu, 10 Aug 2023 12:00:53 -0700 Subject: [PATCH 06/11] removal of unecessary function; fix on the scipt description --- src/compass/utils/correct_cslc_amplitude.py | 13 ++----------- src/compass/utils/mosaic_cslc_backscatter.py | 2 +- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py index 5e683b8c..2bddd0e1 100755 --- a/src/compass/utils/correct_cslc_amplitude.py +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -14,16 +14,6 @@ PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' -def _get_epsg(gdal_raster_path: str): - ''' - Get the EPSG of the input raster - ''' - ds_input = gdal.Open(gdal_raster_path, gdal.GA_ReadOnly) - projection = ds_input.GetProjection() - srs = osr.SpatialReference(wkt=projection) - return int(srs.GetAuthorityCode(None)) - - def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): ''' Load the amplitude from CSLC GDAL dataset. @@ -191,7 +181,8 @@ def get_parser(): Get the parser for CLI ''' parser = argparse.ArgumentParser( - description='Comparison script with burst ID', + description=('Extracts CSLC amplitude with thermal noise and / or ' + 'radiometric correction applied'), formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument('-c', dest='cslc_path', diff --git a/src/compass/utils/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py index 40235701..4b59752a 100755 --- a/src/compass/utils/mosaic_cslc_backscatter.py +++ b/src/compass/utils/mosaic_cslc_backscatter.py @@ -321,7 +321,7 @@ def mosaic_nearest_centroid_mode(cslc_raster_list, local_incidence_angle_list, def compute_mosaic_geotransform_dimension(ds_burst_list): '''Compute the GeoTransform vector that covers all rasters in `ds_burst_list` ''' - + # Initialize list to store extents if len(ds_burst_list) == 0: raise RuntimeError('Empty ds_burst_list was provided') From e6aba15c1112eba880a516f8fb44f669c2420c00 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Thu, 10 Aug 2023 12:10:08 -0700 Subject: [PATCH 07/11] fix on docstring for `correct_cslc_amplitude.py` --- src/compass/utils/correct_cslc_amplitude.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py index 2bddd0e1..ab48c577 100755 --- a/src/compass/utils/correct_cslc_amplitude.py +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -1,7 +1,8 @@ #! python ''' -Script to mosaic the input CSLC rasters +Script to extract amplitude of CSLC, +with thermal noise correction and / or radiometric normalization applied ''' import argparse From c8650102793665cee1a364f958a93ca1a05a8e1b Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Thu, 10 Aug 2023 15:06:19 -0700 Subject: [PATCH 08/11] apply sqrt() to the correction factor --- src/compass/utils/correct_cslc_amplitude.py | 2 +- src/compass/utils/mosaic_cslc_backscatter.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py index ab48c577..945da5c8 100755 --- a/src/compass/utils/correct_cslc_amplitude.py +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -51,7 +51,7 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None print('Applying radiometric mormalization') local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) # Apply radiometric normalization using cotangent(local incidence angle) - correction_factor = (np.sin(local_incidence_angle_arr_rad) + correction_factor = np.sqrt(np.sin(local_incidence_angle_arr_rad) / np.cos(local_incidence_angle_arr_rad)) arr_cslc *= correction_factor diff --git a/src/compass/utils/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py index 4b59752a..5b48869c 100755 --- a/src/compass/utils/mosaic_cslc_backscatter.py +++ b/src/compass/utils/mosaic_cslc_backscatter.py @@ -223,7 +223,7 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None print('Applying radiometric mormalization') local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) # Apply radiometric normalization using cotangent(local incidence angle) - correction_factor = (np.sin(local_incidence_angle_arr_rad) + correction_factor = np.sqrt(np.sin(local_incidence_angle_arr_rad) / np.cos(local_incidence_angle_arr_rad)) arr_cslc *= correction_factor @@ -444,7 +444,7 @@ def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, for i_cslc, cslc_path in enumerate(cslc_path_list): cslc_static_path = cslc_static_path_list[i_cslc] - print('Loading CSLC layer in amplutude, and nouise LUTs: ' + print('Loading CSLC layer in amplutude, and noise LUTs:', f'{i_cslc + 1} / {len(cslc_path_list)}', end='\r') datasets_cslc = get_cslc_gdal_dataset(cslc_path, cslc_static_path, From 284a3f451e1eafd305058b79caeac879c2465aee Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Thu, 10 Aug 2023 15:16:31 -0700 Subject: [PATCH 09/11] revision on correction equaton --- src/compass/utils/correct_cslc_amplitude.py | 3 +-- src/compass/utils/mosaic_cslc_backscatter.py | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py index 945da5c8..f820ae71 100755 --- a/src/compass/utils/correct_cslc_amplitude.py +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -51,8 +51,7 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None print('Applying radiometric mormalization') local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) # Apply radiometric normalization using cotangent(local incidence angle) - correction_factor = np.sqrt(np.sin(local_incidence_angle_arr_rad) - / np.cos(local_incidence_angle_arr_rad)) + correction_factor = np.sqrt(np.tan(local_incidence_angle_arr_rad)) arr_cslc *= correction_factor return arr_cslc diff --git a/src/compass/utils/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py index 5b48869c..2d21ba44 100755 --- a/src/compass/utils/mosaic_cslc_backscatter.py +++ b/src/compass/utils/mosaic_cslc_backscatter.py @@ -223,8 +223,7 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None print('Applying radiometric mormalization') local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) # Apply radiometric normalization using cotangent(local incidence angle) - correction_factor = np.sqrt(np.sin(local_incidence_angle_arr_rad) - / np.cos(local_incidence_angle_arr_rad)) + correction_factor = np.sqrt(np.tan(local_incidence_angle_arr_rad)) arr_cslc *= correction_factor return arr_cslc From aba1db602fca16556c221906f0c946b8bc78b418 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sat, 12 Aug 2023 11:29:33 -0700 Subject: [PATCH 10/11] adjustment for noise LUT in HDF; Removal of mosaicking script --- src/compass/utils/correct_cslc_amplitude.py | 4 +- src/compass/utils/mosaic_cslc_backscatter.py | 571 ------------------- 2 files changed, 2 insertions(+), 573 deletions(-) delete mode 100755 src/compass/utils/mosaic_cslc_backscatter.py diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py index f820ae71..0992a07e 100755 --- a/src/compass/utils/correct_cslc_amplitude.py +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -12,7 +12,7 @@ # Constants for dataset location in HDF5 file PATH_CSLC_LAYER_IN_HDF = '/data' PATH_LOCAL_INCIDENCE_ANGLE = '/data/local_incidence_angle' -PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' +PATH_NOISE_LUT = '/metadata/noise_information/thermal_noise_lut' def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): @@ -97,7 +97,7 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, path_cslc = f'{prefix_netcdf}:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' path_local_incidence = f'NETCDF:{cslc_static_path}:{PATH_LOCAL_INCIDENCE_ANGLE}' - path_noise_lut = f'NETCDF:{cslc_path}:{PATH_NOISE_LAYER_IN_HDF}/thermal_noise_lut' + path_noise_lut = f'NETCDF:{cslc_path}:{PATH_NOISE_LUT}' ds_in = gdal.Open(path_cslc, gdal.GA_ReadOnly) diff --git a/src/compass/utils/mosaic_cslc_backscatter.py b/src/compass/utils/mosaic_cslc_backscatter.py deleted file mode 100755 index 2d21ba44..00000000 --- a/src/compass/utils/mosaic_cslc_backscatter.py +++ /dev/null @@ -1,571 +0,0 @@ -#! python - -''' -Script to mosaic the input CSLC rasters -''' -import argparse -from collections import Counter - -import numpy as np -from osgeo import gdal, osr - -# Constants for dataset location in HDF5 file -PATH_CSLC_LAYER_IN_HDF = '/data' -PATH_LOCAL_INCIDENCE_ANGLE = '/data/local_incidence_angle' -PATH_NOISE_LAYER_IN_HDF = '/metadata/noise_information' - - -def get_most_frequent_epsg(gdal_raster_list: list): - ''' - Get the most frequent EPSG among the rasters in 'gdal_raster_list' - - Parameters - ---------- - gdal_raster_list: list - List of strings for GDAL raster dataset - - Returns - ------- - most_common_epsg: int - Most common EPSG code among the input parameters - - ''' - # Get a list of all EPSG codes - epsg_list = [_get_epsg(raster) for raster in gdal_raster_list] - - # Count the occurrence of each EPSG code - counter = Counter(epsg_list) - - # Return the most common EPSG code - most_common_epsg = counter.most_common(1)[0][0] - - return most_common_epsg - - -def _get_epsg(gdal_raster_path: str): - ''' - Get the EPSG of the input raster - ''' - ds_input = gdal.Open(gdal_raster_path, gdal.GA_ReadOnly) - projection = ds_input.GetProjection() - srs = osr.SpatialReference(wkt=projection) - return int(srs.GetAuthorityCode(None)) - - -def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, epsg_out, dx, dy, snap=0, - noise_correction=True, - radiometric_normalization=True, - resampling_alg='BILINEAR'): - ''' - Get the GDAL dataset for CSLC layer and resample & - reprojecte noise LUT (when user opted in) - - Parameters - ---------- - cslc_path: str - path to the CSLC HDF5 file - cslc_static_path: str - path to the CSLC static layer HDF5 file - pol: str - Polarization - epsg_out: int - EPSG of the outputs - dx, dy: float - x / y spacings of the outputs - snap: float - Snapping value when the input rasters needs to be resampled - noise_correction: bool - Flag to turn on/off noise correction - radiometric_normalization:bool - Flag to turn on/ott radiometric normalization - resampling_alg: str - Resampling algorithm for gdal.Warp() - - Returns - ------- - (ds_cslc, - ds_incidence_angle, - ds_noise_resampled): tuple - Respectively, GDAL raster dataset for CSLC, - GDAL raster dataset for local incidence angle, and - GDAL raster dataset for noise LUT - (resampled to the same geogrid as the other two) - ''' - # Get the geotransform and dimensions - - - prefix_netcdf = 'DERIVED_SUBDATASET:AMPLITUDE:NETCDF' - - path_cslc = f'{prefix_netcdf}:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' - path_local_incidence = f'NETCDF:{cslc_static_path}:{PATH_LOCAL_INCIDENCE_ANGLE}' - path_noise_lut = f'NETCDF:{cslc_path}:{PATH_NOISE_LAYER_IN_HDF}/thermal_noise_lut' - - ds_in = gdal.Open(path_cslc, gdal.GA_ReadOnly) - epsg_cslc = _get_epsg(path_cslc) - - ds_noise_lut = (gdal.Open(path_noise_lut, gdal.GA_ReadOnly) if - noise_correction else None) - - ds_local_incidence_angle = (gdal.Open(path_local_incidence, gdal.GA_ReadOnly) if - radiometric_normalization else None) - - # Prepare for the reprojection - # Reproject the CSLC layer, and - # Define the source and target spatial references - gt_in = ds_in.GetGeoTransform() - xsize = ds_in.RasterXSize - ysize = ds_in.RasterYSize - - src_srs = osr.SpatialReference() - src_srs.ImportFromWkt(ds_in.GetProjectionRef()) - dst_srs = osr.SpatialReference() - dst_srs.ImportFromEPSG(epsg_out) - - # decide whether or not to warp CSLC and static layer - flag_warp_cslc_and_static = ((epsg_cslc != epsg_out) - or (gt_in[0] != dx) - or (abs(gt_in[3]) != abs(dy))) - - # Define a coordinate transformation - transform = osr.CoordinateTransformation(src_srs, dst_srs) - - # Define the corner points - corners = [(0, 0), (0, ysize), (xsize, ysize), (xsize, 0)] - - # Convert the corner points in `epsg_out` - transformed_corners = [] - - # Compute the extet of the resampled raster - if epsg_cslc == epsg_out: - extent_reprojected=[gt_in[0], - gt_in[3] + gt_in[5] * ysize, - gt_in[0] + gt_in[1] * xsize, - gt_in[3]] - else: - for corner in corners: - x_from, y_from = (gt_in[0] + corner[0] * gt_in[1], - gt_in[3] + corner[1] * gt_in[5]) - x_transformed, y_transformed, _ = transform.TransformPoint(x_from, y_from) - transformed_corners.append([x_transformed, y_transformed]) - transformed_corners = np.array(transformed_corners) - - if snap != 0: - transformed_corners = np.round(transformed_corners / snap) * snap - - # Extract the extent of the reprojected raster - extent_reprojected = [np.min(transformed_corners[:, 0]), - np.min(transformed_corners[:, 1]), - np.max(transformed_corners[:, 0]), - np.max(transformed_corners[:, 1])] - - # Put together the warp options - warp_options = gdal.WarpOptions(format='MEM', - resampleAlg=resampling_alg, - xRes=dx, - yRes=abs(dy), - outputBounds=extent_reprojected, - dstSRS=dst_srs) - # Resample the noise LUT - ds_noise_resampled = (gdal.Warp('', ds_noise_lut, options=warp_options) if - ds_noise_lut else None) - - if not flag_warp_cslc_and_static: - # Do not perform reprojection when - # EPSG of the burst is the same as that of the mosaic AND - # the spacing of the CSLC is the same as the mosaic's. - # return the GDAL Dataset objects as they are, except for noise LUT - return (ds_in, ds_local_incidence_angle, ds_noise_resampled) - - ds_incidence_angle_reprojected = \ - (gdal.Warp('', ds_local_incidence_angle, options=warp_options) - if ds_local_incidence_angle else None) - ds_cslc_reprojected = gdal.Warp('', ds_in, options=warp_options) - - # Return the warped dataset - return (ds_cslc_reprojected, - ds_incidence_angle_reprojected, - ds_noise_resampled) - - -def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None): - ''' - Load the amplitude from CSLC GDAL dataset. - Apply noise correction and/or radiometric normalization when - the data are provided. - - Parameters - ---------- - ds_cslc_amp: osgeo.gdal.Dataset - GDAL Raster dataset for CSLC amplitude - ds_local_incidence_angle: osgeo.gdal.Dataset - GDAL Raster dataset for local incidence angle - ds_noise_lut: osgeo.gdal.Dataset - GDAL Raster dataset noise LUT - - Returns - ------- - arr_cslc: np.ndarray - CSLC amplitude - ''' - - # Load CSLC amplutude into array - arr_cslc = ds_cslc_amp.ReadAsArray() - - if ds_noise_lut is not None: - print('Applying noise removal') - # Apply noise correction - - arr_cslc = arr_cslc ** 2 - ds_noise_lut.ReadAsArray() - arr_cslc[arr_cslc<0.0] = 0.0 - arr_cslc = np.sqrt(arr_cslc) - - if ds_local_incidence_angle is not None: - print('Applying radiometric mormalization') - local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) - # Apply radiometric normalization using cotangent(local incidence angle) - correction_factor = np.sqrt(np.tan(local_incidence_angle_arr_rad)) - arr_cslc *= correction_factor - - return arr_cslc - - -def find_upperleft_pixel_index(geotransform_burst, geotransform_mosaic): - ''' - Find the upperleft corner of the burst image in mosaic image, in pixels - - Parameters - ---------- - geotransform_burst: tuple - Geotransform parameter for burst - geotransform_mosaic: tuple - Geotransform parameter for mosaic - - Returns - ------- - upperleft_y_px, upperleft_x_px: int - relative location of upper-left corner of the burst in mosaic - ''' - xmin_burst = geotransform_burst[0] - xmin_mosaic = geotransform_mosaic[0] - - ymax_burst = geotransform_burst[3] - ymax_mosaic = geotransform_mosaic[3] - - upperleft_x_px = int((xmin_burst - xmin_mosaic) / geotransform_mosaic[1]) - upperleft_y_px = int((ymax_burst - ymax_mosaic) / geotransform_mosaic[5]) - - return (upperleft_y_px, upperleft_x_px) - - -def mosaic_list_order_mode(cslc_raster_list, - local_incidence_angle_list, - noise_lut_raster_list, - geotransform_mosaic, - shape_mosaic): - ''' - Perform the mosaicking in "list order mode" - i.e. The pixels from the earlier raster gets replaced by the following rasters - - Parameters - ---------- - cslc_raster_list: list - List of GDAL raster dataset for CSLC - local_incidence_angle_list: list - List of GDAL raster dataset for local incidence angle - noise_lut_raster_list: list - List of GDAL raster dataset for noise LUT - geotransform_mosaic: tuple - Geotramsform parameter for the output mosaic - shape_mosaic: tuple - Shape of the output mosaic array - - Returns - ------- - array_mosaic: np.ndarray - Output mosaic as numpy array - ''' - - array_mosaic = np.zeros(shape_mosaic, dtype=np.float32) / 0 - - for i_raster, cslc_raster in enumerate(cslc_raster_list): - print(f'Processing: {i_raster + 1} of {len(cslc_raster_list)}') - noise_lut_raster = noise_lut_raster_list[i_raster] - local_incidence_angle_raster = local_incidence_angle_list[i_raster] - uly, ulx = find_upperleft_pixel_index(cslc_raster.GetGeoTransform(), - geotransform_mosaic) - amplitude_burst = load_amplitude(cslc_raster, - local_incidence_angle_raster, - noise_lut_raster) - length_burst, width_burst = amplitude_burst.shape - - subset_array_mosaic = array_mosaic[uly : uly + length_burst, - ulx : ulx + width_burst] - - mask_valid = ~np.isnan(amplitude_burst) - subset_array_mosaic[mask_valid] = amplitude_burst[mask_valid] - - return array_mosaic - - -def mosaic_nearest_centroid_mode(cslc_raster_list, local_incidence_angle_list, - noise_lut_raster_list, - geotransform_mosaic, shape_mosaic): - ''' - Perform the mosaicking in "nearest centroid" mode - - ''' - raise NotImplementedError('This mosaic function was not implemented yet') - - -def compute_mosaic_geotransform_dimension(ds_burst_list): - '''Compute the GeoTransform vector that covers all rasters in `ds_burst_list` - ''' - - # Initialize list to store extents - if len(ds_burst_list) == 0: - raise RuntimeError('Empty ds_burst_list was provided') - - extents = [] - for ds_burst in ds_burst_list: - gt_burst = ds_burst.GetGeoTransform() - # Compute extent from GeoTransform - xmin = gt_burst[0] - xmax = gt_burst[0] + gt_burst[1] * ds_burst.RasterXSize - ymin = gt_burst[3] + gt_burst[5] * ds_burst.RasterYSize - ymax = gt_burst[3] - extents.append([xmin, xmax, ymin, ymax]) - - extents = np.array(extents) - - # Find common extent - xmin = np.min(extents[:, 0]) - xmax = np.max(extents[:, 1]) - ymin = np.min(extents[:, 2]) - ymax = np.max(extents[:, 3]) - - # Create and return GeoTransform - gt_burst = ds_burst_list[0].GetGeoTransform() - gt_mosaic = (xmin, gt_burst[1], gt_burst[2], ymax, gt_burst[4], gt_burst[5]) - - width_mosaic = int((xmax - xmin)/ gt_burst[1] + 0.5) - height_mosaic = int((ymin - ymax)/ gt_burst[5] + 0.5) - return gt_mosaic, (width_mosaic, height_mosaic) - - -def save_mosaicked_array(mosaic_arr, geotransform_mosaic, epsg_mosaic, - mosaic_filename): - ''' - Save mosaicked array into GDAL Raster - - Parameters - ---------- - mosaic_arr: np.ndarray - Mosaicked image as numpy array - geotransform_mosaic: tuple - Geotransform parameter for mosaic - epsg_mosaic: - EPSG for mosaic as projection parameter - mosaic_filename: - File name of the mosaic raster to save - ''' - - # Create the projection from the EPSG code - srs_mosaic = osr.SpatialReference() - srs_mosaic.ImportFromEPSG(epsg_mosaic) - projection_mosaic = srs_mosaic.ExportToWkt() - - length_mosaic, width_mosaic = mosaic_arr.shape - driver = gdal.GetDriverByName('GTiff') - ds_out = driver.Create(mosaic_filename, - width_mosaic, length_mosaic, 1, gdal.GDT_Float32, - options=['COMPRESS=LZW', 'BIGTIFF=YES']) - - ds_out.SetGeoTransform(geotransform_mosaic) - ds_out.SetProjection(projection_mosaic) - - ds_out.GetRasterBand(1).WriteArray(mosaic_arr) - - ds_out = None - - -def run(cslc_path_list, cslc_static_path_list, pol, mosaic_path, - dx_mosaic=None, dy_mosaic=None, snap_meters=30, - epsg_mosaic=None, mosaic_mode='list_order', - apply_noise_correcion=True, apply_radiometric_normalization=True, - resampling_alg='BILINEAR'): - ''' - Workflow of the CSLC backscatter mosaic - - Parameters - ---------- - cslc_path_list - pol: str - Polarization to mosaic - mosaic_path: - Path to the output mosaic - dx_mosaic: float, Optional - Spacing of the mosaic in x direction - dy_mosaic: float, Optional - Spacing of the mosaic in y direction - snap_meters: - Snapping in meters - epsg_mosaic: int, optional - EPSG of the output mosaic. - If None, then the EPSG will be determined based on the input bursts - mosaic_mode: - Mosaic mode - apply_noise_correcion: bool - Flag whether or not to apply thermal noise correction - resampling_alg: str - Resampling algorithm, if necessary - ''' - - # Determine the EPSG of the mosaic when it is not provided - if epsg_mosaic is None: - print('EPSG was not provided. ' - 'Finding most common projection among the input rasters.') - epsg_mosaic = get_most_frequent_epsg( - [f'NETCDF:{cslc_path}:{PATH_CSLC_LAYER_IN_HDF}/{pol}' - for cslc_path in cslc_path_list]) - print(f'EPSG of the mosaic will be: {epsg_mosaic}') - - if dx_mosaic is None or dy_mosaic is None: - ds_1st_cslc = gdal.Open(f'NETCDF:{cslc_path_list[0]}:{PATH_CSLC_LAYER_IN_HDF}/{pol}') - gt_1st_cslc = ds_1st_cslc.GetGeoTransform() - dx_mosaic = gt_1st_cslc[1] if dx_mosaic is None else dx_mosaic - dy_mosaic = gt_1st_cslc[5] if dy_mosaic is None else dy_mosaic - - # Get the list of gdal raster dataset - num_raster_in = len(cslc_path_list) - cslc_amp_ds_list = [None] * num_raster_in - noise_lut_ds_list = [None] * num_raster_in - local_incidence_angle_ds_list = [None] * num_raster_in - - for i_cslc, cslc_path in enumerate(cslc_path_list): - cslc_static_path = cslc_static_path_list[i_cslc] - print('Loading CSLC layer in amplutude, and noise LUTs:', - f'{i_cslc + 1} / {len(cslc_path_list)}', - end='\r') - datasets_cslc = get_cslc_gdal_dataset(cslc_path, cslc_static_path, - pol, epsg_mosaic, - dx_mosaic, dy_mosaic, snap_meters, - apply_noise_correcion, - apply_radiometric_normalization, - resampling_alg) - (cslc_amp_ds_list[i_cslc], - local_incidence_angle_ds_list[i_cslc], - noise_lut_ds_list[i_cslc]) = datasets_cslc - - print(' ') - # Get the size of mosaic - gt_mosaic, (width_mosaic, height_mosaic) = \ - compute_mosaic_geotransform_dimension(cslc_amp_ds_list) - - # Compute the array for the mosaicked image - mosaic_functions = { - 'list_order': mosaic_list_order_mode, - 'nearest': mosaic_nearest_centroid_mode - } - - if mosaic_mode not in mosaic_functions: - raise NotImplementedError(f'Mosaicking was not implemented for mode: {mosaic_mode}') - - mosaic_arr = mosaic_functions[mosaic_mode]( - cslc_amp_ds_list, local_incidence_angle_ds_list, noise_lut_ds_list, gt_mosaic, - (height_mosaic, width_mosaic)) - - # Write out the mosaicked image - print(f'Writing mosaic to: {mosaic_path}') - save_mosaicked_array(mosaic_arr, gt_mosaic, epsg_mosaic, mosaic_path) - - - -def get_parser(): - ''' - Get the parser for CLI - ''' - parser = argparse.ArgumentParser( - description='Comparison script with burst ID', - formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('-c', - dest='cslc_list_file', - type=str, - default='', - help='cslc list file') - - parser.add_argument('-s', - dest='cslc_static_list_file', - type=str, - default=[], - help='static list file') - - parser.add_argument('-o', - dest='mosaic_out_path', - type=str, - default=None, - help='path to the burst DB') - - parser.add_argument('-p', - dest='pol', - type=str, - default='VV', - help='Polarization') - - parser.add_argument('-te', - dest='target_spacing', - type=float, - nargs=2, - default=[None, None], - help=('spacing of the mosaic. If not provided, ' - 'the spacing of the first input CSLC will be used.'), - metavar=('NUM1', 'NUM2')) - - parser.add_argument('--snap', - dest='snap', - type=float, - default=30.0, - help='Snapping value of the mosaic grid.') - - parser.add_argument('--mode', - dest='mosaic_mode', - type=str, - default='list_order', - help='Mosaic mode') - - parser.add_argument('--noise_correction_off', - dest='apply_noise_correction', - default=True, - action='store_false', - help='Turn off the noise correction') - - parser.add_argument('--radiometric_normalization_off', - dest='apply_radiometric_normalization', - default=True, - action='store_false', - help='Turn off the radiometric normalization') - - - return parser - - -def main(): - ''' - Entrypoint of the fiunction - ''' - parser = get_parser() - args = parser.parse_args() - - with (open(args.cslc_list_file, 'r') as fin_cslc, - open(args.cslc_static_list_file, 'r') as fin_cslc_static): - list_cslc = fin_cslc.read().rstrip('\n').split('\n') - list_cslc_static = fin_cslc_static.read().rstrip('\n').split('\n') - - run(list_cslc, list_cslc_static, args.pol, args.mosaic_out_path, - dx_mosaic=args.target_spacing[0], dy_mosaic=args.target_spacing[1], - snap_meters=args.snap, mosaic_mode=args.mosaic_mode, - apply_noise_correcion=args.apply_noise_correction, - apply_radiometric_normalization=args.apply_radiometric_normalization) - - -if __name__=='__main__': - main() From 7f8f2b7226a0aa469b8769e1c84c4547e32584d2 Mon Sep 17 00:00:00 2001 From: Seongsu Jeong Date: Sat, 12 Aug 2023 11:38:05 -0700 Subject: [PATCH 11/11] code cleanup and docstrig revision --- src/compass/utils/correct_cslc_amplitude.py | 31 +++++++++++---------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/compass/utils/correct_cslc_amplitude.py b/src/compass/utils/correct_cslc_amplitude.py index 0992a07e..c0371d40 100755 --- a/src/compass/utils/correct_cslc_amplitude.py +++ b/src/compass/utils/correct_cslc_amplitude.py @@ -44,13 +44,14 @@ def load_amplitude(ds_cslc_amp, ds_local_incidence_angle=None, ds_noise_lut=None # Apply noise correction arr_cslc = arr_cslc ** 2 - ds_noise_lut.ReadAsArray() - arr_cslc[arr_cslc<0.0] = 0.0 + arr_cslc[arr_cslc < 0.0] = 0.0 arr_cslc = np.sqrt(arr_cslc) if ds_local_incidence_angle is not None: print('Applying radiometric mormalization') - local_incidence_angle_arr_rad = np.deg2rad(ds_local_incidence_angle.ReadAsArray()) - # Apply radiometric normalization using cotangent(local incidence angle) + local_incidence_angle_arr_rad = \ + np.deg2rad(ds_local_incidence_angle.ReadAsArray()) + # Apply radiometric normalization correction_factor = np.sqrt(np.tan(local_incidence_angle_arr_rad)) arr_cslc *= correction_factor @@ -75,7 +76,7 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, Polarization noise_correction: bool Flag to turn on/off noise correction - radiometric_normalization:bool + radiometric_normalization: bool Flag to turn on/ott radiometric normalization resampling_alg: str Resampling algorithm for gdal.Warp() @@ -85,10 +86,10 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, (ds_cslc, ds_incidence_angle, ds_noise_resampled): tuple - Respectively, GDAL raster dataset for CSLC, - GDAL raster dataset for local incidence angle, and - GDAL raster dataset for noise LUT - (resampled to the same geogrid as the other two) + Respectively, GDAL raster dataset for CSLC, + GDAL raster dataset for local incidence angle, and + GDAL raster dataset for noise LUT + (resampled to the same geogrid as the other two) ''' # Get the geotransform and dimensions @@ -132,8 +133,8 @@ def get_cslc_gdal_dataset(cslc_path, cslc_static_path, pol, dstSRS=src_srs) # Resample the noise LUT - ds_noise_resampled = (gdal.Warp('', ds_noise_lut, options=warp_options) if - ds_noise_lut else None) + ds_noise_resampled = (gdal.Warp('', ds_noise_lut, options=warp_options) + if ds_noise_lut else None) return (ds_in, ds_local_incidence_angle, @@ -148,13 +149,13 @@ def save_amplitude(amplitude_arr, geotransform_cslc, epsg_cslc, Parameters ---------- amplitude_arr: np.ndarray - Mosaicked image as numpy array + Amplitude image as numpy array geotransform_cslc: tuple Geotransform parameter for mosaic - epsg_cslc: + epsg_cslc: int EPSG for mosaic as projection parameter - amplitude_filename: - File name of the mosaic raster to save + amplitude_filename: str + GEOTIFF file name of the output raster ''' # Create the projection from the EPSG code @@ -225,7 +226,7 @@ def get_parser(): def main(): ''' - Entrypoint of the fiunction + Entrypoint of the script ''' parser = get_parser() args = parser.parse_args()