Problems with extraterrestrial bodies
While the great AROSICS software seems to work fine with non-Earth projections in general (at least in global mode), it drops out when using the "verbose" flag in local mode:
(arosics) [swalter]$ arosics local -min_reliability 0 -v 1 -progress 0 hrscref.tif T01_000823_1723_XN_07S080W.tif 100
Calculating footprint polygon and actual data corner coordinates for reference image...
Bounding box of calculated footprint for reference image:
(-4789570.0, -488063.5, -4755145.0, -417426.0)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Bounding box of calculated footprint for image to be shifted:
(-4789363.784239773, -487857.0600601906, -4755351.284239773, -417625.8100601906)
resolutions: 12.5 6.25
100.00 percent of the image to be shifted is covered by the reference image.
Matching window position (X,Y): -4772315.033682144/-452733.69250876445
Initializing tie points grid...
Equalizing pixel grids and projections of reference and target image...
Calculating tie point grid (4543 points) using 32 CPU cores...
Found 1590 matches.
Performing validity checks...
0 tie points flagged by level 1 filtering (reliability).
78 tie points flagged by level 2 filtering (SSIM).
436 tie points flagged by level 3 filtering (RANSAC)
1076 valid tie points remain after filtering.
Visualizing CoReg points grid...
/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/py_tools_ds/geo/projection.py:167: RuntimeWarning: Could not find a suitable EPSG code for the input WKT string.
warnings.warn('Could not find a suitable EPSG code for the input WKT string.', RuntimeWarning)
ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body
ERROR 6: Cannot find coordinate operations from `PROJCRS["Equirectangular MARS",BASEGEOGCRS["GCS_MARS",DATUM["D_MARS",ELLIPSOID["MARS_localRadius",3396000,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Equidistant Cylindrical (Spherical)",ID["EPSG",1029]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]' to `EPSG:4326'
Traceback (most recent call last):
File "/home/swalter/.conda/envs/arosics/bin/arosics", line 11, in <module>
sys.exit(main())
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/arosics/arosics_cli.py", line 364, in main
parsed_args.func(parsed_args)
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/arosics/arosics_cli.py", line 108, in run_local_coreg
CRL.correct_shifts()
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/arosics/CoReg_local.py", line 810, in correct_shifts
self.calculate_spatial_shifts()
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/arosics/CoReg_local.py", line 484, in calculate_spatial_shifts
self.view_CoRegPoints(figsize=(10, 10))
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/arosics/CoReg_local.py", line 556, in view_CoRegPoints
fig, ax = backgroundIm.show_map(figsize=figsize,
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/geoarray/baseclasses.py", line 1365, in show_map
gA2plot = GeoArray(*self._get_plottable_image(xlim, ylim, band, boundsMap=boundsMap,
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/geoarray/baseclasses.py", line 1183, in _get_plottable_image
image2plot, gt, prj = warp_ndarray(image2plot, self.geotransform, self.projection,
File "/home/swalter/.conda/envs/arosics/lib/python3.8/site-packages/py_tools_ds/geo/raster/reproject.py", line 495, in warp_ndarray
raise Exception('Warping Error: ' + gdal.GetLastErrorMsg())
Exception: Warping Error: Cannot find coordinate operations from `PROJCRS["Equirectangular MARS",BASEGEOGCRS["GCS_MARS",DATUM["D_MARS",ELLIPSOID["MARS_localRadius",3396000,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Equidistant Cylindrical (Spherical)",ID["EPSG",1029]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]' to `EPSG:4326'
The error is well known to me and points to the new way the Proj library is performing reprojections (since version 5) and prohibits the use of unknown celestial bodies. Obviously EPSG:4326 is hard-coded somewhere when trying to transform the projected coordinates to lat/lon.
Here you can find the two used files to track the issue: https://box.fu-berlin.de/f/243958069
Any help is greatly appreciated!
Thanks, Sebastian