Incorrect CRS affects on writers.gdal result

Question

I have a .las file with CRS that is not decoded correctly and looks like following:

At the same time the part of metadata which (as I understand) stands for CRS looks like this:

"comp_spatialreference": "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]]]"

The problem is that when I'm trying to create raster by the following pipeline

pipeline = [
        input_path,
        {
            "filename": raster,
            "gdaldriver": "GTiff",
            "resolution": self.resolution,
            "output_type": "mean",
            "type": "writers.gdal",
            "data_type": "float32",
        }
    ]

the MemoryError: bad allocation occurs. However, the files with correctly decoded CRS that 5 times bigger are processed perfectly. At the same time the pipeline

pipeline = [
        {
            "type": "readers.las",
            "filename": input_path,
            "default_srs": "EPSG:9001"
        },
        {
            "type": "filters.reprojection",
            "out_srs": "EPSG:4326"
        },
        {
            "filename": raster,
            "gdaldriver": "GTiff",
            "resolution": self.resolution,
            "output_type": "mean",
            "type": "writers.gdal",
            "data_type": "float32",
        }
    ]

solves problem in some way - it creates raster that looks like black square.

So, the problems are:

1. How to fix the CRS?
2. How to create raster correctly?

Answer 1

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",
        }
    ]