[GIS] QGIS 3.10 not recognising projection information from raster created using GDAL

coordinate systemgdalosgeoqgis-3raster

Problem

QGIS 3.10 only partly reads projection information from raster file however QGIS 3.4 correctly reads information from same raster file.

I thought I had an issue creating the raster – I may still have but thought I should check my project was correctly set up in QGIS and discovered different behaviours between 3.4, 3.10.6 and 3.12.

Is my raster file set up correctly or have I missed something in QGIS?

Have set up QGIS by loading OSTN15 NTV2 but QGIS 3.10 fails to load the raster file as documented below. 3.4 loads same file successfully. Neither versions had the grid files installed. Information on OSTN15 installation below.

Running Raster Layer Unique Values

QGIS 3.4 – Projection: OSGB 1936 / British National Grid (EPSG:27700)

QGIS 3.10 & 3.12 – Raster Layer Unique Values – Projection: Unknown CRS: BOUNDCRS[SOURCECRS[PROJCRS["OSGB 1936 / British Na�

QGIS project steps

Create new project
Set project CRS to epsg:27700 – British National Grid
Load map layer (OS Master map section)
Load raster file
QGIS 3.4 loads immediately but in QGIS 3.10.6.0 I am faced with the window below. Note it states the source CS is unknown 3.4 recognises the source CS.
Multiple operations are possible for converting co-ordinates between these two coordinate reference systems. / Source CRS Unknown CRS
I can select cancel and my image is still drawn in the correct location and correct size. The transformation is not processed as the same error exists when I check the raster layer properties.

Raster Creation Code

src_filename = "000249.png"
dst_filename = "raster249.tif"
fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
dst_ds = driver.Create(dst_filename, xsize=1152, ysize=1152,
                    bands=1, eType=gdal.GDT_Byte)
raster = cv2.imread('000249.png',0)
dst_ds.SetGeoTransform([324643.6636, .017, 0, 673776.9037, 0, -.017])
srs = osr.SpatialReference()
srs.ImportFromEPSG(27700)
dest_wkt = srs.ExportToWkt()
dst_ds.SetProjection(dest_wkt)
dst_ds.GetRasterBand(1).WriteArray(raster)
dst_ds = None

Raster Format

Driver: GTiff/GeoTIFF

Size is 1152 x 1152 x 1

Projection is PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","27700"]]

Origin = (324643.6636, 673776.9037)

Pixel Size = (0.017, -0.017)

Band Type=Byte

Min=0.000, Max=155.000

The raster code was initially obtained from the gdal tutorial but updated from the following stackoverflow questions.

https://gdal.org/tutorials/raster_api_tut.html

Georeferencing raster using GDAL and Python?

ogr: how to use epsg database for RD New projection?

how to set the geotransform using GDAL in python

Loading OSTN15 NTv2 in QGIS

@Ian Turton suggested checking that I have the grid files for epsg:27700 BNG installed which I had not. Here are the steps I followed.

Download zip file from Ordnance Survey
https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/for-developers

I unzipped and copied the files to C:\Program Files\QGIS 3.10\share\proj, there are other .gsb files there. Reading the posts the srs.db file that is released with QGIS should be set up correctly and only the OS files were required.
http://blog.sourcepole.ch/2014/02/18/ntv2-transformations-with-qgis/

https://www.xyht.com/gnsslocation-tech/using-ostn15-qgis

CRS project settings

Best Answer

If I'm not mistaken, it's an issue with the way QGIS (since 3.10) report CRS info taken from GDAL 3.0.4. This issue on GitHub explains everything. What you can do is use this python script written by rouault to update the metadata of your GeoTIFF (it removes the TOWGS84 string that cause the problem). The bug has been patched in GDAL 3.1.0 and when QGIS will start using this version instead of GDAL 3.0.4, your original file should be recognized correctly.

I had the same problem with other CRS (EPSG 2945 to 2952) and using the script did the trick.

Related Solutions

[GIS] CRS transform from EPSG 4326 (WGS 84) ? to EPSG 27700 (OSGB 1936 / British National Grid)

Solved, thanks.

Needed to being ST_Transform instead of ST_SetSRID (see last comment / link)

ALTER TABLE public.roads_sco ALTER COLUMN geom TYPE geometry(MULTILINESTRING, 27700) USING ST_Transform(ST_SetSRID(geom,4326),27700)

Incorrect CRS affects on writers.gdal result

EPSG:9001 may be IGS97 geocentric CRS, or metre Unit. In the WKT definition is used in the nodes of length units, so it isn't the code for the CRS of the dataset.

If we replace the wrongly decoded characters with the ones we think they correspond, and we save the definition in a text file, for example WKT.txt:

COMPD_CS["GCS_Пулково_1942_ГОСТ_32453-2017",PROJCS["GCS_Пулково_1942_ГОСТ_32453-2017",GEOGCS["unknown",DATUM["unnamed",SPHEROID["Krassowsky 1940",6378245,298.300000376014,AUTHORITY["EPSG","7024"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",51],PARAMETER["scale_factor",1],PARAMETER["false_easting",9500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]],VERT_CS["unknown",VERT_DATUM["unknown",2005],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Up",UP]]]

And try to identify it with PROJ:

$ projinfo --identify @WKT.txt
PROJ.4 string:
+proj=tmerc +lat_0=0 +lon_0=51 +k=1 +x_0=9500000 +y_0=0 +a=6378245 +rf=298.300000376014 +units=m +vunits=m +no_defs +type=crs

WKT2:2019 string:
COMPOUNDCRS["GCS_Пулково_1942_ГОСТ_32453-2017",
    PROJCRS["GCS_Пулково_1942_ГОСТ_32453-2017",
        BASEGEOGCRS["unknown",
            DATUM["unnamed",
                ELLIPSOID["Krassowsky 1940",6378245,298.300000376014,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",7024]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]]],
        CONVERSION["unnamed",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",51,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",1,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",9500000,
                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]]]],
    VERTCRS["unknown",
        VDATUM["unknown"],
        CS[vertical,1],
            AXIS["up",up,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]]]

Identification match count: 0

We can see at the end of the output that no match is found for the definition.

Likewise, PROJ converts to the PROJ string and WKT2 formats by default (the input is in WKT1 format).

Therefore, any of these strings can be sent to PDAL so that it assigns the CRS to the dataset when its own definition isn't recognized as valid (or directly does not have it):

pipeline = [
        {
            "type": "readers.las",
            "filename": input_path,
            "default_srs": "+proj=tmerc +lat_0=0 +lon_0=51 +k=1 +x_0=9500000 +y_0=0 +a=6378245 +rf=298.300000376014 +units=m +vunits=m +no_defs +type=crs"
        },
        {
            "type": "filters.reprojection",
            "out_srs": "EPSG:4326"
        },
        {
            "filename": raster,
            "gdaldriver": "GTiff",
            "resolution": self.resolution,
            "output_type": "mean",
            "type": "writers.gdal",
            "data_type": "float32",
        }
    ]