Using geocoded SLC input products to mintpy

[scottstanie]

@taliboliver and I have been working on getting a prep_ script for dolphin and sweets. I'm checking out what metadata attributes would be required/used by mintpy here: https://mintpy.readthedocs.io/en/latest/api/attributes/#required_attributes

• We have X/Y_STEP/START/UNIT for geocoded inputs, but it also saysRANGE_PIXEL_SIZE is required. Our starting point is an X/Y_PIXEL_SIZE for GSLCs
• Similar question about ALOOKS/RLOOKS, which are really XLOOKS/YLOOKS for us

(@taliboliver feel free to add any others that have come up for you)

[taliboliver]

Additionally, the stacks from both dolphin and sweets are in UTM, currently MintPy only recognizes lat/lon besides raster coordinates. What would be the best approach to integrate functionality such as atmospheric corrections or GPS comparisons among others?

[yunjunz]

I am excited to hear you are working on integrating GSLC with mintpy!

RANGE_PIXEL_SIZE is required. Our starting point is an X/Y_PIXEL_SIZE for GSLCs. Similar question about ALOOKS/RLOOKS, which are really XLOOKS/YLOOKS for us.

The A/RLOOKS is required, and the RANGE/AZIMUTH_PIXEL_SIZE is required for radar-coordinate products only. I would still suggest having them in the metadata for completeness. And yes, they are X/Y_PIXEL_SIZE and X/YLOOKS for you. Regarding the attributes name, I would suggest keeping using RANGE/AZIMUTH_PIXEL_SIZE and A/RLOOKS, because 1) range / azimuth is already used indistinguishably with x / y in mintpy when applicable; 2) range / azimuth is more accurate than x / y for airborne data, due to its flexible flight direction.

The ASF HyP3 products are in UTM coordinates, and have been working with mintpy for a while now, except for a few places as marked in https://github.com/insarlab/MintPy/milestone/6. You may find prep_hyp3.py useful as a reference. Besides the usually required metadata, UTM_ZONE and EPSG are required for UTM products, which are extracted from the GDAL VRT file currently for hyp3.

For other steps down the workflow, we could use Y_UNIT == meter with UTM_ZONE or EPSG metadata to identify if a file is in UTM, then apply the changes needed.

1. For InSAR-native steps in mintpy, we should really adjust the code to handle both cases (lat/lon and UTM) in a uniform way as much as possible. This is mostly done while adding support for hyp3 I believe. A few remaining changes could be the template option names, such as mintpy.reference.lalo --> mintpy.reference.en to be more generic.

2. For other functionalities, which mintpy is connected to, e.g. pyaps, gacos, google earth, having a suite of functions to convert 0/1/2D matrix from UTM to WGS84 and vice versa would be needed. Then one could call them in the existing function whenever it's needed. The 0/1/2D conversion from UTM to WGS84 is in utils.utils0.to_latlon(). I would suggest having a utils/utm.py sub-module to handle all the UTM coordinate related conversions.

[yunjunz]

The UTM to/from lat/lon conversion is added as utils.utils0.utm2latlon() and utils.utils0.latlon2utm() in #1052 leveraging the utm module.

[scottstanie]

A couple ongoing notes while putting together a geometry.h5 like https://mintpy.readthedocs.io/en/latest/api/data_structure/#geometry

1. The isce3 convention for heading might be different than isce2?
   For the heading angle, the Compass outputs this (where sweets stitches the bursts into one image) for an ascending track
   
   where looks like this convention for heading
   

I saw that there's the line in prep_hyp3.py

    # ensure negative value for the heading angle
    meta['HEADING'] = float(hyp3_meta['Heading']) % 360. - 360.

which would make this around -260

then there's the line

meta['ORBIT_DIRECTION'] = 'ASCENDING' if abs(meta['HEADING']) < 90 else 'DESCENDING'

which would label it as descending

2. I see theres a latitude (azimuthCoord) and longitude (rangeCoord) there- do you think that's necessary for geocoded products when the whole grid is defined by a X/Y_FIRST, X/Y_STEP?

related to (2), I started on an attempt to add the few change needed to make the HDF5 files CF-compliant so you could load them into xarray/QGIS/gdal. I'm doing this in dolphin here (and probably will soon move these up to a more general module for easier outside use)
https://github.com/opera-adt/dolphin/blob/0ba6b373ce3635322966edc5fbd944c6a62aa246/src/dolphin/workflows/_product.py#L267-L322
basically, for geocoded HDF5 files, we could just add

• x, y 1D datasets with the lon/lon or UTM coordinates of the grid
• use the attach_scale method to make them recognized
• add some metadata to a grid_mapping dummy dataset that GDAL reads

let me know if this path interests you, since we could also do it for the 3D datasets like timeseries.h5

[hfattahi]

The discussion on heading angle is relevant to this issue already open in compass. We are working on providing azimuth angle of the LOS vector which is what we really need. The issue is being worked on upstream in isce3. The other option is to simply provide ENU components.
I assume if we have only ENU components, during loading to mintpy, we have to still compute incidence angle and azimuth angles because that is what Mintpy currently uses. Is that correct?

[yunjunz]

For geometries, mintpy uses the same convention as isce2, which defines HEADING angle from the north clockwise. So yes, it's different from isce3 convention as you describe above.

I assume if we have only ENU components, during loading to mintpy, we have to still compute incidence angle and azimuth angles because that is what Mintpy currently uses. Is that correct?

Right. Maybe because I have not used the unit vector much, the LOS incidence and azimuth angles seem easier, and could be directly used to calculate the unit vector for azimuth offset as well.

I see there's a latitude (azimuthCoord) and longitude (rangeCoord) there- do you think that's necessary for geocoded products when the whole grid is defined by an X/Y_FIRST, X/Y_STEP?

No need. The X/Y_FIRST/STEP/UNIT metadata is enough. For ARIA, the latitude/longitude datasets are not loaded into mintpy, and I have not heard anyone complained.

make the HDF5 files CF-compliant

That sounds really cool @scottstanie! The changes seem small. I am definitely interested. It would make products much more user-friendly.

[scottstanie]

Just an update on how this is going: I got some version working, but wanted to check- Do you have any opinions about where the function/functions should go which take care of this? it's around ~100 lines of code.

The fucntion can be run on the existing timeseries.h5 file to make it recognized by GDAL (as long as you force it to use the NetCDF driver) and make it loadable through xarray using engine="h5netcdf:. We could also add it to the initial creation of the .h5 files too.

$ gdalinfo NETCDF:timeseries_sequential.h5:timeseries

Driver: netCDF/Network Common Data Format
Files: timeseries_sequential.h5
Size is 9636, 6932
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 11N",
    BASEGEOGCRS["WGS 84",
    ....
    ID["EPSG",32611]]
Data axis to CRS axis mapping: 1,2
Origin = (480540.000000000000000,3902670.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  NC_GLOBAL#ALOOKS=1
  .....
  Corner Coordinates:
Upper Left  (  480540.000, 3902670.000) (117d12'50.23"W, 35d16' 1.07"N)
Lower Left  (  480540.000, 3694710.000) (117d12'33.28"W, 33d23'29.41"N)
Upper Right (  769620.000, 3902670.000) (114d 2'14.68"W, 35d13'51.15"N)
Lower Right (  769620.000, 3694710.000) (114d 6' 9.01"W, 33d21'28.26"N)
Center      (  625080.000, 3798690.000) (115d38'25.63"W, 34d19'18.99"N)
Band 1 Block=302x109 Type=Float32, ColorInterp=Undefined
  Metadata:
    DIMENSION_LABELS={time,y,x}
    grid_mapping=spatial_ref
    NETCDF_DIM_time=0
    NETCDF_VARNAME=timeseries
Band 2 Block=302x109 Type=Float32, ColorInterp=Undefined
.... (150 more bands)

ds = xr.open_dataset("timeseries_sequential.h5", engine='h5netcdf')
In [37]: ds = xr.open_dataset("timeseries_sequential.h5", engine='h5netcdf')

In [38]: ds
Out[38]:
<xarray.Dataset>
Dimensions:      (time: 151, y: 6932, x: 9636)
Coordinates:
  * time         (time) datetime64[ns] 2017-02-27 2017-05-10 ... 2020-12-26
  * x            (x) float64 4.806e+05 4.806e+05 ... 7.696e+05 7.696e+05
  * y            (y) float64 3.903e+06 3.903e+06 ... 3.695e+06 3.695e+06
Data variables:
    bperp        (time) float32 ...
    date         (time) |S8 ...
    spatial_ref  int64 ...
    timeseries   (time, y, x) float32 ...
Attributes: (12/26)
    ALOOKS:              1
    AZIMUTH_PIXEL_SIZE:  14.1
    CENTER_LINE_UTC:     49436.654675
    EARTH_RADIUS:        6371000.0
    EPSG:                32611
    FILE_TYPE:           timeseries
    ...                  ...
    X_FIRST:             480540.0
    X_STEP:              30.0
    X_UNIT:              meters
    Y_FIRST:             3902670.0
    Y_STEP:              -30.0
    Y_UNIT:              meters

the HDF5 structure additions looks like this:

$ h5ls -r timeseries.h5
/                        Group
/bperp                   Dataset {151}
/date                    Dataset {151}
/spatial_ref             Dataset {SCALAR}
/time                    Dataset {151}
/timeseries              Dataset {151/Inf, 6932, 9636}
/x                       Dataset {9636}
/y                       Dataset {6932}

• I called the 3rd dimension "time" since "date" was taking with the string representation, and i didn't want to break backward compatibility
• i'm not 100% sure i got all of the NetCDF-specific attributes saved (like _Netcdf4Dimid and _Netcdf4Coordinates that get populated by xarray), but it seems to work for a lot of purposes
• i believe that in order to drag into qgis, you'd need to rename with the extension .nc. i dont know how to force QGIS to not use the HDF5 driver with a .h5 file, even if it's CF-compliant

[yunjunz]

The functionality seems to fit in mintpy.utils.writefile.py, next to the layout_hdf5(), which does most of the HDF5 file creation.

We could somehow call this function within layout_hdf5() to control the engine="h5netcdf: I imagine. And since the change to files is backward compatible, I don't see a problem to enable this in the initial h5 creation yet. One question for this: if we enable the h5netcdf engine, the previous h5 file reading via h5py is not affected, right?

I called the 3rd dimension "time" since "date" was taking with the string representation, and i didn't want to break backward compatibility

Maybe we could use "datetime" instead? Neither "date" nor "time" feels accurate, given that we handle both spaceborne and airborne data.

[scottstanie]

Oh sorry I missed some of these other questions-
The engine="h5netcdf" is something specific to xarray. It tries to use the netcdf library to read .nc files, but you can also specify h5netcdf which works even if you have some partially invalid datasets (like us with the stacks whose depth dimension is (date1, date2).
There's not really a way to specify it in the file; the user opening something with xarray needs to have h5netcdf installed.

Also, i can move the functions in the separate file I made in the PR to be in writefile.py like you suggested