Closed jhkennedy closed 1 month ago
Okay, I think I've got this figured out. The S1 correction workflow only depends on a secondary scene for the dt between the reference and secondary pair. In the HyP3 autoRIFT correction workflow, we only use the reference scene and spoof the dt for the fake secondary scene: https://github.com/ASFHyP3/hyp3-autorift/blob/develop/src/hyp3_autorift/s1_isce2.py#L312-L315
Looking at this pair of S1 scenes, which were acquired almost exactly 12 days apart:
S1A_IW_SLC__1SDV_20170203T162106_20170203T162132_015121_018B9A_1380
S1A_IW_SLC__1SDV_20170215T162105_20170215T162133_015296_019127_6E1E
If I change the m11m12.py
script to divide by 12.0 when loading in the window_rdr_off2vel_*_vec.tif
GeoTIFFs:
offset2vy_1, _, _, _, _ = utils.load_geospatial(offset2vy_path, band=1)
offset2vy_1[offset2vy_1 == nodata] = np.nan
offset2vy_1 /= 12.
offset2vy_2, _, _, _, _ = utils.load_geospatial(offset2vy_path, band=2)
offset2vy_2[offset2vy_2 == nodata] = np.nan
offset2vy_2 /= 12.
offset2vx_1, _, _, _, _ = utils.load_geospatial(offset2vx_path, band=1)
offset2vx_1[offset2vx_1 == nodata] = np.nan
offset2vx_1 /= 12.
offset2vx_2, _, _, _, _ = utils.load_geospatial(offset2vx_path, band=2)
offset2vx_2[offset2vx_2 == nodata] = np.nan
offset2vx_2 /= 12.
offset2vr, _, _, _, _ = utils.load_geospatial(offset2vx_path, band=3)
offset2vr[offset2vr == nodata] = np.nan
offset2vr /= 12.
I can reproduce very close to the correct M11/M12 data values and dr_to_vr_factor
attribute:
short M11(y, x) ;
M11:_FillValue = -32767s ;
M11:dr_to_vr_factor = 70.8575143225911 ;
M11:scale_factor = 0.0002382159f ;
M11:add_offset = 0.009221046f ;
short M12(y, x) ;
M12:_FillValue = -32767s ;
M12:dr_to_vr_factor = 70.8575143225911 ;
M12:scale_factor = 0.0002571183f ;
M12:add_offset = 0.002896907f ;
To correct the M11/M12 variables in an S1 velocity granule, the values computed for M11/M12 and the dr_to_vr_factor via the correction workflow must be divided by the pair separation (in days).
Some of the S1 velocity granules produced for the ITS_LIVE project have incorrect M11/M12 variables (https://github.com/nasa-jpl/its_live_production/issues/18), which either need to be reprocessed or corrected.
Ideally, instead of reprocessing the granules, we'd be able to use the S1 correction workflow to derive M11/M12 using the script below since M11/M12 only depends on the output Geogrid GeoTIFFs.
If I re-process this pair using HyP3 AutoRIFT:
I get M11/M12 variables with these attributes (Note: only showing ones of interest):
and notably, if I run the M11/M12 script below on the Geogrid GeoTIFFs produced when reprocessing the velocity granule, I get the same M11/M12 variables out.
However, if I run the correction workflow using HyP3 AutoRIFT:
I get M11/M12 variables that generally look the same: (conference internet isn't letting me upload the image; will try again at the hotel room tonight)
but are ~1 order of magnitude different:
Comparing the GeoTIFFs produced, we see that only the
window_rdr_off2vel_*_vec.tif
tiffs differ (Golden = Velocity; New = Correction):but also visually look the same (when the colorscale is stretched to +/- 2 sdv for each band).
Script to produce M11/M12 NetCDF file from the S1 Correction workflow output GeoTIFFs:
I've uploaded the script to produce a netCDF File containing the M11/M12 variables derived from output GeoGrid GeoTIFFs to: ``` aws s3 ls --no-sign-request s3://jhk-shenanigans/ITS_LIVE/s1_correction/m11m12/ ``` The files under that prefix are: ``` s3://jhk-shenanigans/ITS_LIVE/s1_correction/m11m12/ ├── correction │ ├── S1A_IW_SLC__1SDV_20170203T162106_20170203T162132_015121_018B9A_1380_conversion_matrices.nc │ ├── window_chip_size_max.tif │ ├── window_chip_size_min.tif │ ├── window_location.tif │ ├── window_offset.tif │ ├── window_rdr_off2vel_x_vec.tif │ ├── window_rdr_off2vel_y_vec.tif │ ├── window_scale_factor.tif │ ├── window_search_range.tif │ └── window_stable_surface_mask.tif ├── m11m12.py └── velocity ├── S1A_IW_SLC__1SDV_20170203T162106_20170203T162132_015121_018B9A_1380_conversion_matrices.nc ├── window_chip_size_max.tif ├── window_chip_size_min.tif ├── window_location.tif ├── window_offset.tif ├── window_rdr_off2vel_x_vec.tif ├── window_rdr_off2vel_y_vec.tif ├── window_scale_factor.tif ├── window_search_range.tif └── window_stable_surface_mask.tif ``` where the `velocity` directory contains GeoGrid GeoTIFFs from the full S1 velocity pair workflow, and the `correction` directory contains GeoGrid GeoTIFFs from the S1 correction workflow. The full script is: ```python from datetime import datetime import numpy as np from autoRIFT import __version__ as version from netCDF4 import Dataset from hyp3_autorift import utils from hyp3_autorift.process import DEFAULT_PARAMETER_FILE scene = 'S1A_IW_SLC__1SDV_20170203T162106_20170203T162132_015121_018B9A_1380' epsg: int = 32635 grid_location_path: str = 'window_location.tif' search_range_path: str = 'window_search_range.tif' offset2vx_path: str = 'window_rdr_off2vel_x_vec.tif' offset2vy_path: str = 'window_rdr_off2vel_y_vec.tif' scale_factor_path: str = 'window_scale_factor.tif' parameter_file: str = DEFAULT_PARAMETER_FILE xGrid, tran, _, srs, nodata = utils.load_geospatial(grid_location_path, band=1) offset2vy_1, _, _, _, _ = utils.load_geospatial(offset2vy_path, band=1) offset2vy_1[offset2vy_1 == nodata] = np.nan offset2vy_2, _, _, _, _ = utils.load_geospatial(offset2vy_path, band=2) offset2vy_2[offset2vy_2 == nodata] = np.nan offset2vx_1, _, _, _, _ = utils.load_geospatial(offset2vx_path, band=1) offset2vx_1[offset2vx_1 == nodata] = np.nan offset2vx_2, _, _, _, _ = utils.load_geospatial(offset2vx_path, band=2) offset2vx_2[offset2vx_2 == nodata] = np.nan offset2vr, _, _, _, _ = utils.load_geospatial(offset2vx_path, band=3) offset2vr[offset2vr == nodata] = np.nan scale_factor_1, _, _, _, _ = utils.load_geospatial(scale_factor_path, band=1) scale_factor_1[scale_factor_1 == nodata] = np.nan dimidY, dimidX = xGrid.shape noDataMask = xGrid == nodata y = np.arange(tran[3], tran[3] + tran[5] * dimidY, tran[5]) x = np.arange(tran[0], tran[0] + tran[1] * dimidX, tran[1]) dr_2_vr_factor = np.median(offset2vr[np.logical_not(np.isnan(offset2vr))]) chunk_lines = np.min([np.ceil(8192 / dimidY) * 128, dimidY]) ChunkSize = [chunk_lines, dimidX] nc_outfile = Dataset(f'{scene}_conversion_matrices.nc', 'w', clobber=True, format='NETCDF4') nc_outfile.setncattr('GDAL_AREA_OR_POINT', 'Area') nc_outfile.setncattr('Conventions', 'CF-1.8') nc_outfile.setncattr('date_created', datetime.now().strftime("%d-%b-%Y %H:%M:%S")) nc_outfile.setncattr('title', 'autoRIFT S1 Corrections') nc_outfile.setncattr('autoRIFT_software_version', version) nc_outfile.setncattr('autoRIFT_parameter_file', parameter_file) nc_outfile.createDimension('x', len(x)) nc_outfile.createDimension('y', len(y)) var = nc_outfile.createVariable('x', np.dtype('float64'), 'x', fill_value=None) var.setncattr('standard_name', 'projection_x_coordinate') var.setncattr('description', 'x coordinate of projection') var.setncattr('units', 'm') var[:] = x var = nc_outfile.createVariable('y', np.dtype('float64'), 'y', fill_value=None) var.setncattr('standard_name', 'projection_y_coordinate') var.setncattr('description', 'y coordinate of projection') var.setncattr('units', 'm') var[:] = y mapping_var_name = 'mapping' var = nc_outfile.createVariable(mapping_var_name, 'U1', (), fill_value=None) if srs.GetAttrValue('PROJECTION') == 'Polar_Stereographic': var.setncattr('grid_mapping_name', 'polar_stereographic') var.setncattr('straight_vertical_longitude_from_pole', srs.GetProjParm('central_meridian')) var.setncattr('false_easting', srs.GetProjParm('false_easting')) var.setncattr('false_northing', srs.GetProjParm('false_northing')) var.setncattr('latitude_of_projection_origin', np.sign(srs.GetProjParm('latitude_of_origin')) * 90.0) var.setncattr('latitude_of_origin', srs.GetProjParm('latitude_of_origin')) var.setncattr('semi_major_axis', float(srs.GetAttrValue('GEOGCS|SPHEROID', 1))) var.setncattr('scale_factor_at_projection_origin', 1) var.setncattr('inverse_flattening', float(srs.GetAttrValue('GEOGCS|SPHEROID', 2))) var.setncattr('spatial_ref', srs.ExportToWkt()) var.setncattr('crs_wkt', srs.ExportToWkt()) var.setncattr('proj4text', srs.ExportToProj4()) var.setncattr('spatial_epsg', epsg) var.setncattr('GeoTransform', ' '.join(str(x) for x in tran)) elif srs.GetAttrValue('PROJECTION') == 'Transverse_Mercator': var.setncattr('grid_mapping_name', 'universal_transverse_mercator') zone = epsg - np.floor(epsg / 100) * 100 var.setncattr('utm_zone_number', zone) var.setncattr('longitude_of_central_meridian', srs.GetProjParm('central_meridian')) var.setncattr('false_easting', srs.GetProjParm('false_easting')) var.setncattr('false_northing', srs.GetProjParm('false_northing')) var.setncattr('latitude_of_projection_origin', srs.GetProjParm('latitude_of_origin')) var.setncattr('semi_major_axis', float(srs.GetAttrValue('GEOGCS|SPHEROID', 1))) var.setncattr('scale_factor_at_central_meridian', srs.GetProjParm('scale_factor')) var.setncattr('inverse_flattening', float(srs.GetAttrValue('GEOGCS|SPHEROID', 2))) var.setncattr('spatial_ref', srs.ExportToWkt()) var.setncattr('crs_wkt', srs.ExportToWkt()) var.setncattr('proj4text', srs.ExportToProj4()) var.setncattr('spatial_epsg', epsg) var.setncattr('GeoTransform', ' '.join(str(x) for x in tran)) else: raise Exception(f'Projection {srs.GetAttrValue("PROJECTION")} not recognized for this program') NoDataValue = -32767 var = nc_outfile.createVariable('M11', np.dtype('int16'), ('y', 'x'), fill_value=NoDataValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) var.setncattr('standard_name', 'conversion_matrix_element_11') var.setncattr( 'description', 'conversion matrix element (1st row, 1st column) that can be multiplied with vx to give range pixel ' 'displacement dr (see Eq. A18 in https://www.mdpi.com/2072-4292/13/4/749)' ) var.setncattr('units', 'pixel/(meter/year)') var.setncattr('grid_mapping', mapping_var_name) var.setncattr('dr_to_vr_factor', dr_2_vr_factor) var.setncattr('dr_to_vr_factor_description', 'multiplicative factor that converts slant range ' 'pixel displacement dr to slant range velocity vr') M11 = offset2vy_2 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) / scale_factor_1 x1 = np.nanmin(M11[:]) x2 = np.nanmax(M11[:]) y1 = -50 y2 = 50 C = [(y2 - y1) / (x2 - x1), y1 - x1 * (y2 - y1) / (x2 - x1)] var.setncattr('scale_factor', np.float32(1 / C[0])) var.setncattr('add_offset', np.float32(-C[1] / C[0])) M11[noDataMask] = NoDataValue * np.float32(1 / C[0]) + np.float32(-C[1] / C[0]) var[:] = M11 var = nc_outfile.createVariable('M12', np.dtype('int16'), ('y', 'x'), fill_value=NoDataValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) var.setncattr('standard_name', 'conversion_matrix_element_12') var.setncattr( 'description', 'conversion matrix element (1st row, 2nd column) that can be multiplied with vy to give range pixel ' 'displacement dr (see Eq. A18 in https://www.mdpi.com/2072-4292/13/4/749)' ) var.setncattr('units', 'pixel/(meter/year)') var.setncattr('grid_mapping', mapping_var_name) var.setncattr('dr_to_vr_factor', dr_2_vr_factor) var.setncattr('dr_to_vr_factor_description', 'multiplicative factor that converts slant range pixel displacement dr to slant range velocity vr') M12 = -offset2vx_2 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) / scale_factor_1 x1 = np.nanmin(M12[:]) x2 = np.nanmax(M12[:]) y1 = -50 y2 = 50 C = [(y2 - y1) / (x2 - x1), y1 - x1 * (y2 - y1) / (x2 - x1)] var.setncattr('scale_factor', np.float32(1 / C[0])) var.setncattr('add_offset', np.float32(-C[1] / C[0])) M12[noDataMask] = NoDataValue * np.float32(1 / C[0]) + np.float32(-C[1] / C[0]) var[:] = M12 ```