Invalid CRS after writing in QgsRasterCalculator

[kadyb]

What is the bug or the crash?

The documentation in QgsRasterCalculator indicates that the CRS will be taken from the first entry in rasterEntries:

outputExtent – output extent. CRS for output is taken from first entry in rasterEntries.

However, in the following example after loading a new raster, the CRS is not valid. Can you confirm that this is a bug?

Steps to reproduce the issue

The sample data is available for download here: https://github.com/kadyb/adg2024/blob/main/dane/DEM.tif

from qgis.core import QgsRasterLayer, QgsCoordinateTransformContext
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

raster = QgsRasterLayer("DEM.tif", "DEM")
print(raster.isValid())
#> True
print(raster.crs())
#> <QgsCoordinateReferenceSystem: EPSG:2180>

a = QgsRasterCalculatorEntry()
a.raster = raster
a.bandNumber = 1
a.ref = "DEM@1"

calculator = QgsRasterCalculator(
    formulaString = "'DEM@1' * 1",
    outputFile = "output.tif",
    outputFormat = "GTiff",
    outputExtent = raster.extent(),
    nOutputColumns = raster.width(),
    nOutputRows = raster.height(),
    rasterEntries = [a],
    transformContext = QgsCoordinateTransformContext()
)

calculator.processCalculation()

test = QgsRasterLayer("output.tif", "test")
print(test.isValid())
#> True
print(test.crs())
#> <QgsCoordinateReferenceSystem: invalid>

Versions

This applies to versions 3.40.0 and 3.34.6 on Windows 10.

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

[kannes]

Is there a sidecar file with an invalid CRS? Try a random name to make 100% sure.

What does gdalinfo say about the output file?

[kadyb]

Is there a sidecar file with an invalid CRS? Try a random name to make 100% sure.

I only have one file in the directory, but anyway, changing the name didn't help.

What does gdalinfo say about the output file?

I included the information below, but it looks like the CRS is just not written to the output file.

Input file

C:\Program Files\QGIS>gdalinfo "C:\Users\Krzysztof\DEM.tif"
Driver: GTiff/GeoTIFF
Files: C:\Users\Krzysztof\DEM.tif
Size is 533, 608
Coordinate System is:
PROJCRS["ETRF2000-PL / CS92",
    BASEGEOGCRS["ETRF2000-PL",
        DATUM["ETRF2000 Poland",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",9702]],
    CONVERSION["Poland CS92",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",19,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9993,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",-5300000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (x)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (medium and small scale)."],
        AREA["Poland - onshore and offshore."],
        BBOX[49,14.14,55.93,24.15]],
    ID["EPSG",2180]]
Data axis to CRS axis mapping: 2,1
Origin = (253698.331199999985984,657570.755200000014156)
Pixel Size = (499.737146529080746,-499.730910526315824)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  253698.331,  657570.755) ( 15d15'54.69"E, 53d43'28.58"N)
Lower Left  (  253698.331,  353734.362) ( 15d29'18.88"E, 50d59'51.47"N)
Upper Right (  520058.230,  657570.755) ( 19d18'16.23"E, 53d46'56.78"N)
Lower Right (  520058.230,  353734.362) ( 19d17'10.49"E, 51d 3' 0.25"N)
Center      (  386878.281,  505652.558) ( 17d20'11.76"E, 52d24'18.29"N)
Band 1 Block=533x3 Type=Float32, ColorInterp=Gray
  NoData Value=9999

Output file

C:\Program Files\QGIS>gdalinfo "C:\Users\Krzysztof\output.tif"
Driver: GTiff/GeoTIFF
Files: C:\Users\Krzysztof\output.tif
Size is 533, 608
Origin = (253698.331199999985984,657570.755200000014156)
Pixel Size = (499.737146529080746,-499.730910526315824)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  253698.331,  657570.755)
Lower Left  (  253698.331,  353734.362)
Upper Right (  520058.230,  657570.755)
Lower Right (  520058.230,  353734.362)
Center      (  386878.281,  505652.558)
Band 1 Block=533x3 Type=Float32, ColorInterp=Gray
  Min=13.000 Max=373.785
  Minimum=13.000, Maximum=373.785, Mean=116.336, StdDev=43.161
  NoData Value=-3.4028235e+38
  Metadata:
    STATISTICS_APPROXIMATE=YES
    STATISTICS_MAXIMUM=373.78491210938
    STATISTICS_MEAN=116.33588384848
    STATISTICS_MINIMUM=13
    STATISTICS_STDDEV=43.160723362679
    STATISTICS_VALID_PERCENT=91.2

[agiudiceandrea]

It seems to me that using the QgsRasterCalculator construction variant without the transformContext parameter (although deprecated since QGIS 3.8)

calculator = QgsRasterCalculator(
    formulaString = "'DEM@1' * 1",
    outputFile = "output.tif",
    outputFormat = "GTiff",
    outputExtent = raster.extent(),
    nOutputColumns = raster.width(),
    nOutputRows = raster.height(),
    rasterEntries = [a]
)

the output raster has a valid CRS

print(test.isValid())
#> True
print(test.crs())
#> <QgsCoordinateReferenceSystem: EPSG:2180>