The discrepancy in latitude and longitude when displayed in QGIS

[kyutecdaichi]

I am currently conducting InSAR analysis using ISCE2-MintPy with Sentinel data.
However, when I convert the geo_temporalCoherence.h5 or geo_velocity.h5 files generated by MintPy to GeoTIFF format (.tif) and display them in QGIS, there is a slight misalignment in the images.
I am unsure about the cause or how to resolve this issue. If anyone knows the reason or has a solution, please let me know. I would greatly appreciate your help.
Thank you!

• Operating system: Ubuntu 18.04.6 LTS
• Python environment: conda
• MintPy version: 1.5.3
• InSAR processor/product: isce2

[welcome]

👋 Thanks for opening your first issue here! Please filled out the template with as much details as possible. We appreciate that you took the time to contribute!
Make sure you read our contributing guidelines.

[codeautopilot]

Potential solution

The plan to solve the bug involves verifying and correcting the geotransform parameters and the coordinate reference system (CRS) used when converting the MintPy-generated HDF5 files to GeoTIFF format. The discrepancy in latitude and longitude alignment in QGIS is likely due to incorrect geotransform settings or an inappropriate EPSG code being applied. By ensuring that the metadata accurately reflects the spatial properties of the data and that the correct CRS is used, we can resolve the misalignment issue.

What is causing this bug?

The bug is likely caused by incorrect geotransform parameters or an inappropriate EPSG code being applied during the conversion process. The geotransform parameters (X_FIRST, X_STEP, Y_FIRST, Y_STEP) define the spatial resolution and origin of the data. If these values are incorrect, the spatial alignment in QGIS will be off. Additionally, if the EPSG code used does not match the actual CRS of the data, it can lead to misalignment. The default to EPSG 4326 (WGS84) may not be suitable for the dataset, especially if it was collected in a different projection.

Code

To address the issue, we need to ensure that the geotransform parameters and CRS are correctly set. Here is a code snippet that demonstrates how to verify and set these parameters:

from osgeo import gdal, osr

def convert_to_geotiff(input_h5, output_tif, metadata):
    # Open the input HDF5 file
    dataset = gdal.Open(input_h5, gdal.GA_ReadOnly)
    
    # Extract geotransform parameters from metadata
    x_first = metadata.get('X_FIRST', None)
    x_step = metadata.get('X_STEP', None)
    y_first = metadata.get('Y_FIRST', None)
    y_step = metadata.get('Y_STEP', None)
    
    if None in [x_first, x_step, y_first, y_step]:
        raise ValueError("Geotransform parameters are missing from metadata.")
    
    # Set geotransform
    geotransform = (x_first, x_step, 0, y_first, 0, y_step)
    
    # Create output GeoTIFF
    driver = gdal.GetDriverByName('GTiff')
    out_dataset = driver.CreateCopy(output_tif, dataset, 0)
    out_dataset.SetGeoTransform(geotransform)
    
    # Set the correct CRS
    srs = osr.SpatialReference()
    epsg_code = metadata.get('EPSG', None)
    if epsg_code:
        srs.ImportFromEPSG(int(epsg_code))
    else:
        raise ValueError("EPSG code is missing from metadata.")
    
    out_dataset.SetProjection(srs.ExportToWkt())
    out_dataset.FlushCache()

    print(f"GeoTIFF saved to {output_tif} with geotransform and CRS set correctly.")

# Example usage
metadata = {
    'X_FIRST': 100.0,
    'X_STEP': 0.0001,
    'Y_FIRST': 0.0,
    'Y_STEP': -0.0001,
    'EPSG': '32633'  # Example EPSG code
}
convert_to_geotiff('geo_temporalCoherence.h5', 'output.tif', metadata)

How to replicate the bug

1. Use MintPy to generate geo_temporalCoherence.h5 or geo_velocity.h5 files from Sentinel data.
2. Convert these HDF5 files to GeoTIFF format using a script that does not correctly set the geotransform parameters or CRS.
3. Open the resulting GeoTIFF files in QGIS.
4. Observe the misalignment in latitude and longitude when compared to other correctly aligned datasets.

By following these steps, you should be able to replicate the misalignment issue in QGIS.

Click here to create a Pull Request with the proposed solution

Files used for this task:

Changes on src/mintpy/save_gdal.py

Based on the provided code in src/mintpy/save_gdal.py, here are some observations and potential areas that could be causing the discrepancy in latitude and longitude when the GeoTIFF files are displayed in QGIS:

1. Geotransform Parameters:

   • The geotransform is set using X_FIRST, X_STEP, Y_FIRST, and Y_STEP from the metadata. If these values are incorrect or not accurately representing the spatial resolution and origin, it could lead to misalignment in QGIS.
   • Ensure that these metadata values are correctly calculated and represent the actual spatial properties of the data.

2. Coordinate Reference System (CRS):

   • The code attempts to set the EPSG code for the spatial reference system using either EPSG or UTM_ZONE from the metadata. If neither is present, it defaults to EPSG 4326 (WGS84).
   • If the default EPSG 4326 is not appropriate for the dataset, this could cause misalignment. Verify that the correct EPSG code is being used.

3. Data Type Conversion:

   • The code converts boolean data to uint8 because GDAL does not support boolean types. While this is unlikely to affect geolocation, ensure that this conversion does not inadvertently affect the data representation.

4. Projection and Transform Setting:

   • The projection is set using the EPSG code, and the transform is applied to the raster. If there is any mismatch between the data's actual projection and the one being set, it could lead to discrepancies.
   • Double-check that the projection and transform settings are consistent with the data's intended geospatial reference.

5. Metadata Integrity:

   • Ensure that the metadata being read and used for geotransform and projection is accurate and complete. Any missing or incorrect metadata could lead to issues in spatial alignment.

In summary, the potential cause of the discrepancy could be related to incorrect geotransform parameters, an inappropriate EPSG code, or metadata issues. It would be beneficial to verify these aspects in the metadata of the input files and ensure they are correctly applied during the conversion process.

Disclaimer: This comment was entirely generated using AI. Be aware that the information provided may be incorrect.

Current plan usage: 0.00%

Have feedback or need help?
Discord
Documentation
[email protected]

[yunjunz]

@kyutecdaichi could you check the alignment of the KMZ file in Google Earth? If there is misalignment there as well, then the problem is in the data itself, due to the inaccurate geolocation, probably from the reference image used in coregistration while running ISCE2/topsStack. If there is no misalignment there, then they may be a problem in the code itself, while saving the HDF5 to GeoTIFF.

[kyutecdaichi]

@yunjunz Thank you for your response.
When I checked the KMZ file on Google Maps, I observed that the same misalignment appeared in the KMZ file as well.
If I reanalyze the data with ISCE2, could you please advise on any points I should pay attention to, such as parameters for stackSentinel or other related settings? Your help would be much appreciated.
Thank you.

[yunjunz]

Use a reference date in 2019 for coregistration, to avoid potential ionospheric displacement.

How big is this misalignment you observed?

[kyutecdaichi]

When I displayed the geo_velocity.kmz on Google Maps, I noticed that the roads were slightly misaligned.
By the way, the analysis period was from 2015 to 2023 (excluding December, January, and February of each year).

[yunjunz]

Is it possible that you are using ALOS-2, processed with isce2/stripmapStack, and forgot to run stackStripMap.py with --zero optioin?

[kyutecdaichi]

I am using Sentinel-1 data, not ALOS-2 data. Therefore, I am using the isce2/topsStack's stackSentinel.py for processing.

[yunjunz]

Interesting. Could you post your basic data info here?

• path number
• bounding box in SNWE
• stackSentinel.py command used for isce2 processing?

[kyutecdaichi]

I am not very familiar with 'path number,' but when checking on ASF data search, should I use Path 46 as shown in the photo?

・path number
　46
・bounding box in SNWE
　'43.031439 43.279942 141.276939 141.652175'
・stackSentinel.py command
　-n '2 3' -c 5
・dem.py command
　dem.py -a stitch -b 42 43 141 142 -r -s 1 -c