2

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

enter image description here

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?
4
  • Hi! It is a Transverse Mercator centered at 51 degrees of longitude with 9500000 false easting, from Kassowsky 1940 ellipsoid, compounded with an unknown vertical reference. EPSG:9001 is the code for meters unit of measure, it isn't the CRS of the dataset. the absence of vertical reference could be a problem, but we can find the WKT definition of that system to be assigned to the original dataset. Aug 16 at 1:41
  • 1
    "default_srs": "+proj=tmerc +lat_0=0 +lon_0=51 +k=1 +x_0=9500000 +y_0=0 +a=6378245 +rf=298.3 +units=m +vunits=m +no_defs +type=crs" might solve the issue. Identification process with PROJ doesn't return any match. The string has WKT1 format, but you can send the PROJ string instead. Please, check if it works. Aug 16 at 2:40
  • Yes, this worked, but there are files without any metadata except name. In this case crs recognition is impossible, right? Aug 16 at 4:06
  • 1
    A default CRS can be assigned in the same way. If the name does not allow us to infer if it is the correct CRS, then we will not know if it is. Aug 16 at 13:18
2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.