search

yumorishita / LiCSBAS

LiCSBAS: InSAR time series analysis package using LiCSAR products
https://doi.org/10.3390/RS12030424
GNU General Public License v3.0
228 stars 109 forks source link

transform of interferograms and DEM are different #135

Closed Fanchengyan closed 2 years ago

Fanchengyan commented 2 years ago

I found that the upper-left corner of interferograms and DEM are half a pixel away from each other. For example, when opening an interferogram file like '20141024_20150128.geo.unw' and the DEM file in this frame 'xxxD_xxxxx_xxxxxx.geo.hgt.tif' in QGIS, there would be a half-pixel shift between each other.

And in the script LiCSBAS02_ml_prep.py, the upper-left corner of interferograms would shift half-pixel before writing this information into EQA.dem_par.

I want to make sure that whether this half-pixel offset is due to the processing errors or the fact that there is indeed a half-pixel offset between interferograms and DEM?

yumorishita commented 2 years ago

Could you show me an example of the half-pixel shift? I cannot find the shift in my data. There should not be the half-pixel offset between interferograms and DEM.

The shift in LiCSBAS02_ml_prep.py is to convert the pixel registration (GeoTIFF) to the grid registration (EQA.dem_par) (e.g., see http://geophysics.eas.gatech.edu/classes/Intro_GMT/gmt_www/gmt/doc/html/GMT_Docs/node151.html). The coordinates in GeoTIFF and EQA.dem_par have the half-pixel difference, but they mean exactly the same locations.

Fanchengyan commented 2 years ago

I tried to load some interferograms and DEMs in different frames into QGIS, they all have the problem of half-pixel shift ( see the screenshot below). I am thinking about the reason for this. Is this due to the processing error that mistakingly treats "PixelIsArea" Raster Space as "PixelIsPoint" Raster Space? http://web.archive.org/web/20160326194152/http://remotesensing.org/geotiff/spec/geotiff2.5.html#2.5.2

For the tiff files, I found it is in the "PixelIsArea" Raster Space. And the coordinate value of west/north (G[0] an G[3] for gdal Geotransform : https://gdal.org/tutorials/geotransforms_tut.html) should be referred to the upper-left corner of the first pixel(upper-left pixel)

If it is due to the above case, I want to know which raster (interferogram or DEM) has the correct/original Geotransform?

截图录屏_选择区域_20211124133747 截图录屏_选择区域_20211124134027 截图录屏_选择区域_20211124135045

yumorishita commented 2 years ago

Could you show me the return of gdalinfo for both unw.tif and hgt.tif files in a frame?

Fanchengyan commented 2 years ago

for DEM:

gdalinfo 106D_05049_131313.geo.hgt.tif -json

{
  "description":"106D_05049_131313.geo.hgt.tif",
  "driverShortName":"GTiff",
  "driverLongName":"GeoTIFF",
  "files":[
    "106D_05049_131313.geo.hgt.tif",
    "106D_05049_131313.geo.hgt.tif.aux.xml"
  ],
  "size":[
    3476,
    2674
  ],
  "coordinateSystem":{
    "wkt":"GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]",
    "dataAxisToSRSAxisMapping":[
      2,
      1
    ]
  },
  "geoTransform":[
    97.703887939453125,
    0.0009999992325902,
    0.0,
    40.8277778625488281,
    0.0,
    -0.0009999992325902
  ],
  "metadata":{
    "":{
      "AREA_OR_POINT":"Point",
      "TIFFTAG_DATETIME":"2016:12:14 12:33:01",
      "TIFFTAG_SOFTWARE":"Created with GAMMA Software www.gamma-rs.ch"
    },
    "IMAGE_STRUCTURE":{
      "COMPRESSION":"LZW",
      "INTERLEAVE":"BAND"
    }
  },
  "cornerCoordinates":{
    "upperLeft":[
      97.7038879,
      40.8277779
    ],
    "lowerLeft":[
      97.7038879,
      38.1537799
    ],
    "lowerRight":[
      101.1798853,
      38.1537799
    ],
    "upperRight":[
      101.1798853,
      40.8277779
    ],
    "center":[
      99.4418866,
      39.4907789
    ]
  },
  "wgs84Extent":{
    "type":"Polygon",
    "coordinates":[
      [
        [
          97.7038879,
          40.8277779
        ],
        [
          97.7038879,
          38.1537799
        ],
        [
          101.1798853,
          38.1537799
        ],
        [
          101.1798853,
          40.8277779
        ],
        [
          97.7038879,
          40.8277779
        ]
      ]
    ]
  },
  "bands":[
    {
      "band":1,
      "block":[
        3476,
        1
      ],
      "type":"Float32",
      "colorInterpretation":"Gray",
      "min":1013.923,
      "max":5575.115,
      "minimum":1013.923,
      "maximum":5575.115,
      "mean":2263.014,
      "stdDev":1194.097,
      "noDataValue":0.0,
      "metadata":{
        "":{
          "STATISTICS_APPROXIMATE":"YES",
          "STATISTICS_MAXIMUM":"5575.1147460938",
          "STATISTICS_MEAN":"2263.0139799297",
          "STATISTICS_MINIMUM":"1013.9225463867",
          "STATISTICS_STDDEV":"1194.0970169161",
          "STATISTICS_VALID_PERCENT":"86.27"
        }
      }
    }
  ]
}

for interferogram:

gdalinfo 20141024_20141117.geo.unw.tif -json
{
  "description":"20141024_20141117.geo.unw.tif",
  "driverShortName":"GTiff",
  "driverLongName":"GeoTIFF",
  "files":[
    "20141024_20141117.geo.unw.tif"
  ],
  "size":[
    3476,
    2674
  ],
  "coordinateSystem":{
    "wkt":"GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]",
    "dataAxisToSRSAxisMapping":[
      2,
      1
    ]
  },
  "geoTransform":[
    97.703388900394998,
    0.00099999921,
    0.0,
    40.8282777996049973,
    0.0,
    -0.00099999921
  ],
  "metadata":{
    "":{
      "AREA_OR_POINT":"Point",
      "TIFFTAG_DATETIME":"2019:11:10 02:04:32",
      "TIFFTAG_SOFTWARE":"Created with GAMMA Software www.gamma-rs.ch data2geotiff v2.7"
    },
    "IMAGE_STRUCTURE":{
      "COMPRESSION":"LZW",
      "INTERLEAVE":"BAND"
    }
  },
  "cornerCoordinates":{
    "upperLeft":[
      97.7033889,
      40.8282778
    ],
    "lowerLeft":[
      97.7033889,
      38.1542799
    ],
    "lowerRight":[
      101.1793862,
      38.1542799
    ],
    "upperRight":[
      101.1793862,
      40.8282778
    ],
    "center":[
      99.4413875,
      39.4912789
    ]
  },
  "wgs84Extent":{
    "type":"Polygon",
    "coordinates":[
      [
        [
          97.7033889,
          40.8282778
        ],
        [
          97.7033889,
          38.1542799
        ],
        [
          101.1793862,
          38.1542799
        ],
        [
          101.1793862,
          40.8282778
        ],
        [
          97.7033889,
          40.8282778
        ]
      ]
    ]
  },
  "bands":[
    {
      "band":1,
      "block":[
        3476,
        1
      ],
      "type":"Float32",
      "colorInterpretation":"Gray",
      "noDataValue":0.0,
      "metadata":{
      }
    }
  ]
}

It's obvious that geoTransform of those two are different

yumorishita commented 2 years ago

I have confirmed that the two files you showed have half-pixel shift. I realised that the shift does not exist in all frames and the shift seems to depend on the creation time of the files. Old (before 2019?) files have consistent coordinates, and new ones as well.

This shift is caused in the LiCSAR processing, not in LiCSBAS. I cannot identify which coordinate is true. I have contacted the people in charge and they will post something.

For your information, the coordinates in LiCSBAS are extracted from an unw GeoTIFF file.

espiritocz commented 2 years ago

Hi.

Many thanks for finding about such issue! To help us find effective solution, can you please confirm that this resampled file fits correctly? https://gws-access.jasmin.ac.uk/public/nceo_geohazards/LiCSAR_products/106/106D_05049_131313/interferograms/20141024_20141117/20141024_20141117.geo.unw.warped.tif

cheers!

Fanchengyan commented 2 years ago

Hi @espiritocz , from what I can see, this resampled file has no problem, and the two images overlap perfectly after your modification.

So, does this mean that the geotransform of DEM is correct? Could I directly assign the geotransform of DEM to Interferogram metric?

espiritocz commented 2 years ago

thank you for the confirmation. the geotiffs must be resampled first.

please redownload the input geotiff files, they are now fixed towards the DEM. We will investigate this issue further to find the optimal solution.

all best!

yumorishita commented 2 years ago

@espiritocz Has this issue been solved?

espiritocz commented 2 years ago

hi Yu,

we solved the issue in LiCSAR but some frames might still be corrupt. please let me know should the error reappear.

On Sun, 31 Jul 2022, 14:42 Yu Morishita, @.***> wrote:

@espiritocz https://github.com/espiritocz Has this issue been solved?

— Reply to this email directly, view it on GitHub https://github.com/yumorishita/LiCSBAS/issues/135#issuecomment-1200417677, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADYYKV3VSJB3WBSMSASZLLVWZYEVANCNFSM5IQD6LQA . You are receiving this because you were mentioned.Message ID: @.***>