Warp Reproject Runs but does nothing

[Zuline]

What is the bug or the crash?

Running MacOS (Apple Silicon M3)

Running Warp Reproject runs OK but the output is unchanged. Running the same command in gdalwarp in terminal results in the desired reprojection.

Steps to reproduce the issue

1. Open Warp Reproject
2. Select a raster layer
3. Enter source and target SRS
4. Specify filepath
5. Run the reprojection

Expected Outcome: File will be reprojected into the target SRS
Actual Outcome: File is created but not reprojected

Versions

QGIS version	3.38.2-Grenoble	QGIS code revision	130c432
Qt version	5.15.2
Python version	3.9.5
GDAL/OGR version	3.3.2
PROJ version	8.1.1
EPSG Registry database version	v10.028 (2021-07-07)
GEOS version	3.9.1-CAPI-1.14.2
SQLite version	3.35.2
PDAL version	2.3.0
PostgreSQL client version	unknown
SpatiaLite version	5.0.1
QWT version	6.1.6
QScintilla2 version	2.11.5
OS version	macOS 14.6
Active Python plugins
landxmlconvertor	0.2.0
QuickOSM	2.2.3
quick_map_services	0.19.34
DEMto3D	3.6
ViewshedAnalysis	1.9
qgis2web	3.22.0
processing_r	4.1.0
qfieldsync	v4.10.1
Qgis2threejs	2.7.3
mmqgis	2021.9.10
mapflow	2.6.2
processing	2.12.99
grassprovider	2.12.99
db_manager	0.1.20
MetaSearch	0.3.6
HCMGIS	24.8.27

QGIS version 3.38.2-Grenoble QGIS code revision [130c432](https://github.com/qgis/QGIS/commit/130c432394e) Qt version 5.15.2 Python version 3.9.5 GDAL/OGR version 3.3.2 PROJ version 8.1.1 EPSG Registry database version v10.028 (2021-07-07) GEOS version 3.9.1-CAPI-1.14.2 SQLite version 3.35.2 PDAL version 2.3.0 PostgreSQL client version unknown SpatiaLite version 5.0.1 QWT version 6.1.6 QScintilla2 version 2.11.5 OS version macOS 14.6

Active Python plugins
landxmlconvertor
0.2.0
QuickOSM
2.2.3
quick_map_services
0.19.34
DEMto3D
3.6
ViewshedAnalysis
1.9
qgis2web
3.22.0
processing_r
4.1.0
qfieldsync
v4.10.1
Qgis2threejs
2.7.3
mmqgis
2021.9.10
mapflow
2.6.2
processing
2.12.99
grassprovider
2.12.99
db_manager
0.1.20
MetaSearch
0.3.6
HCMGIS
24.8.27

Supported QGIS version

• I'm running a supported QGIS version according to the roadmap.

New profile

• I tried with a new QGIS profile

Additional context

No response

[nirvn]

@Zuline , I could not reproduce this issue here, it might be dataset specific? A sample dataset to use alongside your steps to replicate the issue would be needed here.

[Zuline]

@Zuline , I could not reproduce this issue here, it might be dataset specific? A sample dataset to use alongside your steps to replicate the issue would be needed here.

The ortho is too big to attach here. It can be downloaded at this link: https://ln5.sync.com/dl/402dcdf20/zv37byji-8n7k2ecr-43hatcmb-cyixk86p

Regards
Mike

[agiudiceandrea]

I cannot replicate the issue on my Windoww 10 system using QGIS 3.38.2 (GDAL 3.9.2). @Zuline , please specify all the algorithm's parameters used and the full processing log. How do you say the output raster layer is not reprojected? Please also provide the output raster layer.

[Zuline]

I cannot replicate the issue on my Windoww 10 system using QGIS 3.38.2 (GDAL 3.9.2). @Zuline , please specify all the algorithm's parameters used and the full processing log. How do you say the output raster layer is not reprojected? Please also provide the output raster layer.

The input CRS is specified by Warp Reproject as EPSG:32755 from inspection of the ortho. I accept that.
Output CRS required is EPSG:7855 which I specify.
Resampling method is bilinear.

The log is below and the ortho that results is in the same link as above.

The reason I say that it is not reprojected is:

1. The original and the "reprojected" ortho exactly overlay each other - nothing changes, except that the "reprojected" ortho in the information tab now says it's EPSG:7855;
2. I have GCPs with points collected in EPSG:7855. The points when created in QGIS are not coincident on the targets in the original or the "reprojected" ortho - the difference is systematic and around 1.54m
3. When I use gdalwarp at the command line the ortho is clearly reprojected and the known points and the GCP targets are exactly coincident.

`QGIS version: 3.38.2-Grenoble
QGIS code revision: 130c432
Qt version: 5.15.2
Python version: 3.9.5
GDAL version: 3.3.2
GEOS version: 3.9.1-CAPI-1.14.2
PROJ version: Rel. 8.1.1, September 1st, 2021
PDAL version: 2.3.0 (git-version: Release)
Algorithm started at: 2024-09-10T10:33:45
Algorithm 'Warp (reproject)' starting…
Input parameters:
{ 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : '/Users/fred/Library/CloudStorage/Sync/Drone/Mapping/Quarries-1/Output/all/odm_orthophoto/odm_orthophoto.tif', 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : '/Users/fred/Library/CloudStorage/Sync/Drone/Mapping/Quarries-1/reproj-test.tif', 'RESAMPLING' : 1, 'SOURCE_CRS' : QgsCoordinateReferenceSystem('EPSG:32755'), 'TARGET_CRS' : QgsCoordinateReferenceSystem('EPSG:7855'), 'TARGET_EXTENT' : None, 'TARGET_EXTENT_CRS' : None, 'TARGET_RESOLUTION' : None }

GDAL command:
gdalwarp -overwrite -s_srs EPSG:32755 -t_srs EPSG:7855 -r bilinear -of GTiff /Users/fred/Library/CloudStorage/Sync/Drone/Mapping/Quarries-1/Output/all/odm_orthophoto/odm_orthophoto.tif /Users/fred/Library/CloudStorage/Sync/Drone/Mapping/Quarries-1/reproj-test.tif
GDAL command output:
Using band 4 of source image as alpha.
Creating output file that is 12107P x 8234L.
Processing /Users/fred/Library/CloudStorage/Sync/Drone/Mapping/Quarries-1/Output/all/odm_orthophoto/odm_orthophoto.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.
Process completed successfully
Execution completed in 17.58 seconds
Results:
OUTPUT: /Users/fred/Library/CloudStorage/Sync/Drone/Mapping/Quarries-1/reproj-test.tif

Loading resulting layers
Algorithm 'Warp (reproject)' finished`

[agiudiceandrea]

@Zuline, QGIS actually uses gdalwarp (in your case, GDAL/OGR version 3.3.2) to perform the task with the specified parameters and the coordinate conversion is actually performed by the PROJ library (in your case, version 8.1.1) based on the EPSG Registry database version v10.028 (2021-07-07) (in your case) definitions.
What is the GDAL/OGR (and PROJ) versions you use at the command line? Do you use an NTv2 / grid shift file / ... to perform the transformation from GDA2020 to WGS84?
What is the datum transformation set in the project / user profile for the EPSG:32755 / EPSG:7855 transformation?
The EPSG:32755 and EPSG:7855 CRSs have the same definition (apart from the different ellipsoid), so the transformation between the two is almost a null transformation if no Helmert 7 parameters or NTv2 / grid shift is set for the datum transformation.

The provided files odm_orthophoto.tif and reproj-test.tif have different CRSs and extents (you can check using gdalinfo or in the layer properies in QGIS), so it seems to me the processing algorithm "Warp (Reproject)" actually performs the transformation (but probably not the way you expect):

gdalinfo odm_orthophoto.tif
Driver: GTiff/GeoTIFF
Files: odm_orthophoto.tif
Size is 10355, 7346
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 55S",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 55S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",147,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 144┬░E and 150┬░E, southern hemisphere between 80┬░S and equator, onshore and offshore. Australia. Papua New Guinea."],
        BBOX[-80,144,0,150]],
    ID["EPSG",32755]]
Data axis to CRS axis mapping: 1,2
Origin = (324193.596915444941260,5815516.980700185522437)
Pixel Size = (0.019998290184675,-0.019999999656037)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  PREDICTOR=2
Corner Coordinates:
Upper Left  (  324193.597, 5815516.981) (145d 0'12.05"E, 37d47'26.92"S)
Lower Left  (  324193.597, 5815370.061) (145d 0'11.93"E, 37d47'31.68"S)
Upper Right (  324400.679, 5815516.981) (145d 0'20.52"E, 37d47'27.06"S)
Lower Right (  324400.679, 5815370.061) (145d 0'20.39"E, 37d47'31.83"S)
Center      (  324297.138, 5815443.521) (145d 0'16.22"E, 37d47'29.37"S)

gdalinfo reproj-test.tif
Driver: GTiff/GeoTIFF
Files: reproj-test.tif
Size is 12107, 8234
Coordinate System is:
PROJCRS["GDA2020 / MGA zone 55",
    BASEGEOGCRS["GDA2020",
        DATUM["Geocentric Datum of Australia 2020",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",7844]],
    CONVERSION["Map Grid of Australia zone 55",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",147,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Australia - onshore and offshore between 144┬░E and 150┬░E."],
        BBOX[-50.89,144,-9.23,150.01]],
    ID["EPSG",7855]]
Data axis to CRS axis mapping: 1,2
Origin = (324172.128683218965307,5815527.262876833789051)
Pixel Size = (0.019998173459327,-0.019998173459327)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  324172.129, 5815527.263) (145d 0'11.19"E, 37d47'26.57"S)
Lower Left  (  324172.129, 5815362.598) (145d 0'11.04"E, 37d47'31.91"S)
Upper Right (  324414.247, 5815527.263) (145d 0'21.08"E, 37d47'26.74"S)
Lower Right (  324414.247, 5815362.598) (145d 0'20.94"E, 37d47'32.08"S)
Center      (  324293.188, 5815444.930) (145d 0'16.06"E, 37d47'29.33"S)

Maybe @rouault know if there are knows issues.

[Zuline]

@agiudiceandrea thanks for your detailed comments. You are far more across the intimacies of this than I am.

The version at the command line is: GDAL 3.9.2, released 2024/08/13. I don't explicitly use any NtV2 or grid shift files.

The one thing I will say is that it's a fundamental misapprehension to treat this as a null transformation, that is to essentially cast EPSG:32755 and EPSG:7855 as being the same. At some point of equivalence when 7855 was first released that may have been the case or nearly the case. It's not so now. There is now a very substantial difference due to the plate fixed nature of 7855 and the transformation needs to deal with that, as it does at the command line.

[rouault]

The one thing I will say is that it's a fundamental misapprehension to treat this as a null transformation, that is to essentially cast EPSG:32755 and EPSG:7855 as being the same. At some point of equivalence when 7855 was first released that may have been the case or nearly the case. It's not so now. There is now a very substantial difference due to the plate fixed nature of 7855 and the transformation needs to deal with that, as it does at the command line.

That would be more a debate for a discussion on the PROJ mailing list, but QGIS does what PROJ says about that, and PROJ says what the EPSG database says about, and in EPSG, EPSG:32755 is "WGS 84 / UTM zone 55S", based on the generic WGS 84 datum ensemble, not a particular realization of it. And in that case, the transformation between WGS 84 and GDA2020 is a null one. If you worke with a UTM projected CRS based on ITRF2014 (or "WGS 84 (G1762)" instead of WGS 84, then you could get a time-dependent transformation with PROJ.

Demo:

Given test.wkt containing

PROJCRS["ITRF2014 / UTM zone 55S",
    BASEGEOGCRS["ITRF2014",
        DYNAMIC[
            FRAMEEPOCH[2010]],
        DATUM["International Terrestrial Reference Frame 2014",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",9000]],
    CONVERSION["UTM zone 55S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",147,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]]]

$ projinfo -s @test.wkt -t EPSG:7855 -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of UTM zone 55S + Conversion from ITRF2014 (geog2D) to ITRF2014 (geocentric) + ITRF2014 to GDA2020 (1) + Conversion from GDA2020 (geocentric) to GDA2020 (geog2D) + Map Grid of Australia zone 55, 0.03 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.

PROJ string:
+proj=pipeline
  +step +inv +proj=utm +zone=55 +south +ellps=GRS80
  +step +proj=cart +ellps=GRS80
  +step +proj=helmert +x=0 +y=0 +z=0 +rx=0 +ry=0 +rz=0 +s=0 +dx=0 +dy=0 +dz=0
        +drx=0.00150379 +dry=0.00118346 +drz=0.00120716 +ds=0 +t_epoch=2020
        +convention=coordinate_frame
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=utm +zone=55 +south +ellps=GRS80

But to get a shift you need to specify a transformation coordinate epoch:

At 2020.0, as the central epoch of the above transformation is 2020.0, you get no change:

$ echo  500000 5000000 0 2020 | bin/cs2cs "$(cat test.wkt)" EPSG:7855
500000.00	5000000.00 -0.00 2020

At 2024.0, you get a few centimeters of shift:

$ echo  500000 5000000 0 2024 | bin/cs2cs "$(cat test.wkt)" EPSG:7855
499999.95	4999999.78 0.00 2024

But time-dependent coordinate transformations are (I believe, but I may be wrong) not yet available through QGIS UI tools. You may need to use GDAL command line utilities to do that. Cf https://gdal.org/en/latest/user/coordinate_epoch.html#support-in-utilities

[agiudiceandrea]

@rouault, it looks like the issue seems also due to the fact that gdalwarp used through the QGIS processing "Warp (Reproject)" GDAL algoritm on macOS (it uses GDAL 3.3.2 / PROJ 8.1.1) outputs a raster reportedly transformed in a not expected way, while using the same gdalwarp command on a separate command line environment (with GDAL 3.9.2, on the same system, not sure how it was installed and which PROJ version is linked to and which EPSG database and grids uses if any) outputs a raster reportedly transformed in an expected way.

time-dependent coordinate transformations are (I believe, but I may be wrong) not yet available through QGIS UI tools. You may need to use GDAL command line utilities to do that.

The "Warp (Reproject)" processing algorithm allows to specify any additional command line parameter, thus also -s_coord_epoch and -t_coord_epoch if needed.
See also #57719.

It also looks like the reported issue about "Warp Reproject does nothing" is not occurring even on the user system, since the georeferencing of the processing algorithm's output is actually different from the georeferencing of the processing algorithm's input, so Warp (Reproject) actually does something.

Anyway it doesn't look like a QGIS issue.

[rouault]

GDAL algoritm on macOS (it uses GDAL 3.3.2 / PROJ 8.1.1) outputs a raster reportedly transformed in a not expected way,

There are actually several transformations between WGS 84 and GDA2020 : one is a null transformation, and another one assumes that WGS 84 = GDA94, and thus applies a GDA94 <--> GDA2020 datum transformation, which is about 1.8 meter shift. What PROJ selects as the default transformation has evolved over PROJ / EPSG versions. That is the likely cause for different results on different platforms.

Confirmed actually:

• with PROJ master:

$ bin/projinfo -s "WGS 84" -t GDA2020 --spatial-test intersects -o PROJ --single-line
Candidate operations found: 3
-------------------------------------
Operation No. 1:

INVERSE(EPSG):8450, Inverse of GDA2020 to WGS 84 (2), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.

PROJ string:
+proj=noop

-------------------------------------
Operation No. 2:

EPSG:9690, WGS 84 to GDA2020 (3), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 3:

DERIVED_FROM(EPSG):9691, WGS 84 to GDA2020 (4), 3.0 m, Australia - Australian Capital Territory; New South Wales; Northern Territory; Queensland; South Australia; Tasmania; Western Australia; Victoria., at least one grid missing

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=au_icsm_GDA94_GDA2020_conformal_and_distortion.tif +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

Grid au_icsm_GDA94_GDA2020_conformal_and_distortion.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/au_icsm_GDA94_GDA2020_conformal_and_distortion.tif

vs

• with PROJ 8.1.1:

$ ~/proj/install-proj-8.1.1-cmake/bin/projinfo -s "WGS 84" -t GDA2020 --spatial-test intersects -o PROJ --single-line
Candidate operations found: 3
-------------------------------------
Operation No. 1:

EPSG:9690, WGS 84 to GDA2020 (3), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 2:

INVERSE(EPSG):8450, Inverse of GDA2020 to WGS 84 (2), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.

PROJ string:
+proj=noop

-------------------------------------
Operation No. 3:

DERIVED_FROM(EPSG):9691, WGS 84 to GDA2020 (4), 3.0 m, Australia - Australian Capital Territory; New South Wales; Northern Territory; Queensland; South Australia; Tasmania; Western Australia; Victoria., at least one grid missing

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=au_icsm_GDA94_GDA2020_conformal_and_distortion.tif +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

Grid au_icsm_GDA94_GDA2020_conformal_and_distortion.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/au_icsm_GDA94_GDA2020_conformal_and_distortion.tif

As you can notice EPSG publishes all 3 transformations with the same accuracy of 3 m...

[nyalldawson]

I think we should close this. It's ultimately user error, even if the underlying cause is very tricky to understand.