Open kautlenbachs opened 1 year ago
Current geoposition transformation was implemented following gdaltransform source code - https://gdal.org/programs/gdaltransform.html
Currently shelved, this needs a lot of implementation added since DEM, SAR metadata calculations are all in degrees and WGS84.
Transforming WGS84 coordinates to UTM 43N for example.
GDAL transformer:
SNAP calculated ones:
GDAL default transfomer 'OGRCoordinateTransformation' gets invalid results for UTM meter ones. Although this is suggested as per manual - https://gdal.org/tutorials/osr_api_tut.html Instead a 'gdaltransform' utility source code was explored and now 'GDALCreateGenImgProjTransformer2()' is used as per manual - https://web.archive.org/web/20180702062214/http://www.gdal.org:80/warptut.html and https://gdal.org/programs/gdaltransform.html
In SNAP GeoTools package/library is used which is also officially certified by OGC. But already the first conversion gets different results from GDAL transformer. In additions a mathematical transform is applied.
Concluding the above, it would be hard to deliver exactly same conversion results as in SNAP.
There is strategy that we plan to implement. Calculate everything as is in WGS84. And finally warp the raster on GPU with resample operation onto the calculated grid. Grid boundaries/rectangle is calculated using GDAL tranformer thus pixels are resampled appropriately.
But this approach would not produce good results due to the curvature not taken account -> more GCPs needed. Check gdalwarp how to continue with this approach.
More precise dataset bounds can be estimated with "GDALSuggestedWarpOutput2()" - https://gis.stackexchange.com/questions/455055/coordinate-transformation-with-gdal-replicating-geotools-java-package-transform?noredirect=1#comment743074_455055 It has been experimented together with GDAL official source code (apps - gdaltransformer.cpp gdalwarp_lib.cpp) and can replicate:
Upper Left ( 588461.655, 4452379.350) ( 73d45'49.45"W, 50d 4'27.38"S) Lower Left ( 588461.655, 4236227.483) ( 73d42'39.21"W, 52d 1' 3.11"S) Upper Right ( 744706.220, 4452379.350) ( 71d34'58.79"W, 50d 1'49.91"S) Lower Right ( 744706.220, 4236227.483) ( 71d26'14.12"W, 51d58'14.42"S) Center ( 666583.937, 4344303.417) ( 72d37'27.76"W, 51d 1'43.14"S)
which are not that far away from SNAP's calculated boundaries: Corner Coordinates: Upper Left ( 588458.695, 4452379.281) ( 73d45'49.60"W, 50d 4'27.39"S) Lower Left ( 588458.695, 4236220.841) ( 73d42'39.36"W, 52d 1' 3.33"S) Upper Right ( 744693.415, 4452379.281) ( 71d34'59.43"W, 50d 1'49.94"S) Lower Right ( 744693.415, 4236220.841) ( 71d26'14.77"W, 51d58'14.65"S) Center ( 666576.055, 4344300.061) ( 72d37'28.15"W, 51d 1'43.26"S)
Implementing gdalwarp functionality, which would simply "resample" the result from WGS84 to anything else is not reasonable, since it could be invoked as an after processing step and would not use GPU for this anyway. Also interpolating already calculated results (which also results from various interpolation of pixels)
Here is a comparison between snap produced one and gdalwarped
That is bad pixels = 31591
bad pixels = 0.018492%
bad pixels = 184.920162 ppm
avg rel diff = 434.370609%
avg rel diff = 4343706.089542530477047 ppm, cnt = 79734117, pct = 46.672932
median = -648421.049118 ppm
Also could be that wherever coordinate conversion is needed the pipeline is split up. There would be workflow where many kernels are created by steps where WGS84 coordinates are output and then on the host side they are transformed with GDAL and then another kernel is executed and so on and so forth....