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:
- How to fix the CRS?
- How to create raster correctly?
Best Answer
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:
And try to identify it with PROJ:
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):