LAS not compliant to specification after spTransform

[Jean-Romain]

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
las2 = spTransform(las, sp::CRS(SRS_string = "EPSG:26918"))
las_check(las2)
#> 
#>  Checking the data
#>   - Checking coordinates... ✓
#>   - Checking coordinates type... ✓
#>   - Checking coordinates quantization...
#>     ✗ 81589 X coordinates were not stored with a resolution compatible with scale factor 0.01 and offset 214261.661487569
#>     ✗ 81589 Y coordinates were not stored with a resolution compatible with scale factor 0.01 and offset 5021517.24154598
#>   - Checking attributes type... ✓
#>   - Checking ReturnNumber validity... ✓
#>   - Checking NumberOfReturns validity... ✓
#>   - Checking ReturnNumber vs. NumberOfReturns... ✓
#>   - Checking RGB validity... ✓
#>   - Checking absence of NAs... ✓
#>   - Checking duplicated points... ✓
#>   - Checking degenerated ground points... ✓
#>   - Checking attribute population...
#>    ⚠ 'PointSourceID' attribute is not populated.
#>    ⚠ 'EdgeOfFlightline' attribute is not populated.
#>   - Checking gpstime incoherances ✓
#>   - Checking flag attributes... ✓
#>   - Checking user data attribute... ✓
#>  Checking the header
#>   - Checking header completeness... ✓
#>   - Checking scale factor validity... ✓
#>   - Checking point data format ID validity... ✓
#>   - Checking extra bytes attributes validity... ✓
#>   - Checking the bounding box validity...
#>     ✗ Stored resolution of 'Min X' and/or 'Max X' not compatible with 'X offset' and 'X scale factor'
#>     ✗ Stored resolution of 'Min Y' and/or 'Max Y' not compatible with 'Y offset' and 'Y scale factor'
#>   - Checking coordinate reference sytem... ✓
#>  Checking header vs data adequacy
#>   - Checking attributes vs. point format... ✓
#>   - Checking header bbox vs. actual content... ✓
#>   - Checking header number of points vs. actual content... ✓
#>   - Checking header return number vs. actual content... ✓
#>  Checking preprocessing already done 
#>   - Checking ground classification... yes
#>   - Checking normalization... yes
#>   - Checking negative outliers... ✓
#>   - Checking flightline classification... no

[wiesehahn]

If I understand correctly, this issue adresses the problem discussed in https://gis.stackexchange.com/questions/385720/r-sptransform-in-lidr-error-non-quantizable-value and should be fixed since version 3.1.2

with the example data this works (e.g.)

LASfile <- system.file("extdata", "example.laz", package="rlas")
las = readLAS(LASfile)
st_crs(las)$Name
st_bbox(las)
tlas <- sf::st_transform(las, sf::st_crs(26918))
st_crs(tlas)$Name
st_bbox(tlas)

But with own data (lidR version 4.0.3) I still get error message Error: Non quantizable value outside the range of representable values of type 'int'

LASfile <- "own.las"
las <- readLAS(LASfile)
projection(las) <- 4326

# Scale factor X Y Z:       1e-07 1e-07 0.01 
# Offset X Y Z:             0 0 0 
# min X Y Z:                10.32643 51.82899 331.63 
# max X Y Z:                10.64345 51.87805 799.5 

tlas <- sf::st_transform(las, sf::st_crs(25832))

Is this behaviour expected due to possibly errorneous data or a bug on lidR / sf side?

[Jean-Romain]

It should work yet bugs are highly probable. I need a file to reproduce. The four points of the bounding box should be sufficient.

[wiesehahn]

> las_check(las)

 Checking the data
  - Checking coordinates... ✓
  - Checking coordinates type... ✓
  - Checking coordinates range... ✓
  - Checking coordinates quantization... ✓
  - Checking attributes type... ✓
  - Checking ReturnNumber validity...
    ⚠ Invalid data: 904039 points with a return number equal to 0 found.
  - Checking NumberOfReturns validity...
    ⚠ Invalid data: 904039 points with a number of returns equal to 0 found.
  - Checking ReturnNumber vs. NumberOfReturns... ✓
  - Checking RGB validity... ✓
  - Checking absence of NAs... ✓
  - Checking duplicated points... ✓
  - Checking degenerated ground points... skipped
  - Checking attribute population...
    🛈 'PointSourceID' attribute is not populated
    🛈 'ScanDirectionFlag' attribute is not populated
    🛈 'EdgeOfFlightline' attribute is not populated
  - Checking gpstime incoherances
    ✗ 196336 pulses (points with the same gpstime) have points with identical ReturnNumber
  - Checking flag attributes... ✓
  - Checking user data attribute... ✓
 Checking the header
  - Checking header completeness... ✓
  - Checking scale factor validity... ✓
  - Checking point data format ID validity... ✓
  - Checking extra bytes attributes validity... ✓
  - Checking the bounding box validity... ✓
  - Checking coordinate reference system... ✓
 Checking header vs data adequacy
  - Checking attributes vs. point format... ✓
  - Checking header bbox vs. actual content... ✓
  - Checking header number of points vs. actual content... ✓
  - Checking header return number vs. actual content... ✓
 Checking coordinate reference system...
  - Checking if the CRS was understood by R... ✓
 Checking preprocessing already done 
  - Checking ground classification... no
  - Checking normalization... no
  - Checking negative outliers... ✓
  - Checking flightline classification... no
 Checking compression
  - Checking attribute compression...
   -  ScanDirectionFlag is compressed
   -  EdgeOfFlightline is compressed
   -  Classification is compressed
   -  Synthetic_flag is compressed
   -  Keypoint_flag is compressed
   -  Withheld_flag is compressed
   -  ScanAngleRank is compressed
   -  UserData is compressed
   -  PointSourceID is compressed

> projection(las) <- 4979

> st_bbox(las)
    xmin     ymin     xmax     ymax 
10.48973 51.85762 10.51178 51.86423 

> summary(las)
class        : LAS (v1.2 format 1)
memory       : 37.9 Mb 
extent       : 10.48973, 10.51178, 51.85762, 51.86423 (xmin, xmax, ymin, ymax)
coord. ref.  : WGS 84 
area         : 0.95 kunits²
points       : 904 thousand points
density      : 0.96 points/units²
File signature:           LASF 
File source ID:           0 
Global encoding:
 - GPS Time Type: GPS Week Time 
 - Synthetic Return Numbers: no 
 - Well Know Text: CRS is GeoTIFF 
 - Aggregate Model: false 
Project ID - GUID:        00000000-0000-0000-0000-000000000000 
Version:                  1.2
System identifier:         
Generating software:      LidarTools, IDL 8.5.1 
File creation d/y:        194/2022
header size:              227 
Offset to point data:     626 
Num. var. length record:  3 
Point data format:        1 
Point data record length: 28 
Num. of point records:    904039 
Num. of points by return: 0 0 0 0 0 
Scale factor X Y Z:       1e-07 1e-07 0.01 
Offset X Y Z:             0 0 0 
min X Y Z:                10.48973 51.85762 546.33 
max X Y Z:                10.51178 51.86423 1324.13 
Variable Length Records (VLR):
   Variable Length Record 1 of 3 
       Description: GeoKeyDirectoryTag 
       Tags:
          Key 1024 value 2 
          Key 1025 value 0 
          Key 1026 value 0 
          Key 2048 value 4326 
          Key 2049 value 44 
          Key 2052 value 9001 
          Key 2054 value 9102 
          Key 2057 value 0 
          Key 2058 value 1 
          Key 3073 value 50 
          Key 4099 value 9001 
          Key 4096 value 5030 
          Key 4097 value 94 
          Key 3072 value 4979 
   Variable Length Record 2 of 3 
       Description: GeoKeyDoubleParamsTag 
       data:                 6378137 6356752 
   Variable Length Record 3 of 3 
       Description: GeoASCIIParamsTag 
       data:                 LAS file produced by XXXWGS-84LAS file produced by XXXWGS84 Ellipsoid 
Extended Variable Length Records (EVLR):  void

> tlas <- sf::st_transform(las, sf::st_crs(25832))
Error: Non quantizable value outside the range of representable values of type 'int'

Is this sufficient or do you need a file?

[Jean-Romain]

1. I see that your coordinates are recorded with an accuracy of 10-7 degrees
2. You are trying to convert from long lat ranging from 10° to 50° to projected crs. Your point cloud is basically 1km². It is impossible to store a point cloud covering 1km² with an accuracy of 10^-7 (100 nm). You must adjust the scale factor. Maybe try with a centimeter accuracy. This is approximately what 10^-7 degrees corresponds to.

sf::st_transform(las, sf::st_crs(25832), scale = 0.01)

I do agree that this could be detected automatically by lidR

[wiesehahn]

Not sure if I understand you correct. Latitute ranges from 10.48973 to 10.51178 ° and Longitude from 51.85762 to 51.86423 °, so this point cloud covers more or less 1.5 x 1 km.

sf::st_transform(las, sf::st_crs(25832), scale = 0.01) works without error, but what is scale doing? I cant find documentation in the context of st_transform. Does this mean that coordinates are rounded to e.g. 10.490 (xmin)?

[wiesehahn]

Ah ok, so the reported min/max values by e.g. st_bbox are truncated and the values are stored internally with a precision of 10^7 (reported by summary(las) )?
So values would be rounded to e.g. 10.48973 from 10.4897300, right?

[Jean-Romain]

R always truncates the display of numbers when printed.

1/3
#> [1] 0.3333333
options(digits = 12)
1/3
#> [1] 0.333333333333
options(digits = 22)
1/3
#> [1] 0.3333333333333333148296

Also the min/max may not be truncated. Maybe the real min is 10.4897300

[Jean-Romain]

Not sure if I understand you correct. Latitute ranges from 10.48973 to 10.51178 ° and Longitude from 51.85762 to 51.86423 °, so this point cloud covers more or less 1.5 x 1 km.

My bad. I read it wrongly. It was 4 AM here ;-). It does not change anything to my answer except your point cloud does not cover half the world.

sf::st_transform(las, sf::st_crs(25832), scale = 0.01) works without error, but what is scale doing? I cant find documentation in the context of st_transform. Does this mean that coordinates are rounded to e.g. 10.490 (xmin)?

https://gis.stackexchange.com/questions/360247/rescaling-and-reoffsetting-a-point-cloud-with-lidr/360253#360253

Your target coordinates (eg. 60654.23) will be rounded to a centimeter accuracy