Azimuthal Equidistant Projection produces terrain corrected data with a shift

Microwave Toolbox
#1

I was just testing SNAP 5 and tried to do terrain correction directly to the projection that we use in-house. Unfortunately after Terrain correction there is a substantial shift in the output product. As can be seen in the screenshot of QGIS with the Google Satellite Layer below. I’ve used the two graphs AO_Cal.xml -> TC_Equi7.xml in succession to get to the output file. Before viewing the data in QGIS I did a gdal_translate Sigma0_VH.img Sigma0_VH.tif using gdal 2.1.2.

AO_Cal.xml (2.2 KB)
TC_Equi7.xml (3.6 KB)

blob.jpg1861×969 730 KB

Custom azimuthal equidistant projection wrong? (Reprojection Module)
#2

As far as I can tell these projection issues could be solved by using a more recent geotools version.

#3

Updating the geotools version is being worked on but, I believe it’s not straight forward.

#4

I’m doing this right now. I try to use latest GeoTools version. As Luis said it is not straight forward. For example new issues are introduced, at some places we need to change the implementation to adapt for changes and at other places it behaves differently and we need to find out why. Maybe it is a bug in GeoTools maybe, maybe an issue on our side, etc.
But the next release will probably use the latest GeoTools version.

2 Likes
#5

Thanks. The problem is hopefully with Geotools. See https://osgeo-org.atlassian.net/browse/GEOT-2310 for example.

If you have a test implementation on any Github branch I could test it for our use case.

#6

If you really like to test the branch your are very welcome. The branch is named geotools-16.1_mp. It is the same on all senbox repositories. For example https://github.com/senbox-org/snap-engine/tree/geoTools-16.1_mp

There is still a known issue with the rendering of Pins and GCPs.

#7

Unfortunately I still see the same problem, strange.

Some additional observations:

  • I tried to export the BEAM/DIMAP product into a GeoTIFF with SNAP to make sure that gdal is not introducing the shift. But I got the following error

    SEVERE [org.openide.util.RequestProcessor]: Error in RequestProcessor org.netbeans.modules.progress.ui.RunOffEDTImpl$3
java.lang.IllegalArgumentException: Unable to map projectionAzimuthal Equidistant

  • Maybe the build requirements should be updated, first I tried building with maven 3.0.5 which lead to failures maven 3.3.9 worked.

#8

Thanks for the hint regarding the build requirements. I was already on 3.3.9. So I didn’t noticed it.
On master it still works with 3.0.5.

Regarding the projection issue.
Probably the shift is introduced by us. GeoTools has support only since September last year for AzimuthalEquidistant projection.
We have this implementation since August 2014.
I’ve deactivated our implementation. Now the one from GeoTools should be used. Maybe this solves it.

#9

I’ve compared the values in our unit level test with the ones provided by cs2cs. They are equal.
Target CRS: +proj=aeqd +a=6370997.0+b=6370997.0 +lat_0=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
source cords:
0.0 ; 0.0
-180.0 ; 90.0
-173.0 ; -45.5
-75.0 ; -35.0
96.5; -21.0

Output is:
0.0;0.0
-7.50445202206e-26;10007538.6856
-1772868.47834;-14803428.1626
-7000586.40014;-5074782.38583
9964308.84278;-3849686.13682

#10

That seems to solve the problem. The big shift is now gone. :grinning:

Did you also test yours on WGS84 or only on a sphere?

You also have to know that proj.4 did not have the accurate implementation of the projection algorithm until a colleague of mine added it in version 4.9.1. But that only produces an error of a few meters and not the kind of shift we saw.

From our internal tests these are correct values calculated by pyproj.

input coordinate pairs

+proj=longlat +datum=WGS84 +no_defs

-31.627336;30.306273
-14.589038;-43.880131
79.423313;-35.261658
23.456413;10.457987

output coordinate pairs

+proj=aeqd +ellps=WGS84 +lat_0=8.5 +lon_0=21.5 +x_0=5621452.01998 +y_0=5990638.42298

484363.5635332661;9118738.557172634
2306098.2802775027;51.98795392457396
11357485.715251137;755174.4295843598
5835685.017902998;6207782.488239396
#11

Yes, our projection was only defined for the sphere case. This is where the shift came from. But with the new GeoTools impl. it works also for the ellipsoid.

#12

Hi, Thanks again for the fix.

From what I could see on Github this fix was not yet released.
Do you have any idea when it will be included in an official SNAP release?

#13

There is no fixed date yet, but it should be out in July.
There are still some issues to be solved with the new GeoTools version.
The new version fixes old issues but also introduces new ones.