Error when running pipeline on VPC at another location

[wiesehahn]

When creating a Virtual Point Cloud file at another location than the referenced data and then trying to run a pipeline on this VPC I get an error.

E.g.:

f <- system.file("extdata", "MixedConifer.las", package="lasR")

create_vpc <- function(input_files, output_file){
  if (!file.exists(output_file)){    
    pipeline = lasR::write_vpc(output_file)
    lasR::exec(pipeline, on = input_files)
  } else {
    print("VPC already exists!")
  }
}

create_vpc(input_files = f, output_file = "C:/Users/jwiesehahn/Downloads/temp/vpc_lasr.vpc")

exec(reader_las(), on = "C:/Users/jwiesehahn/Downloads/temp/vpc_lasr.vpc")

produces

ERROR: ERROR: cannot open file 'C:\Users\jwiesehahn\Downloads\temp'
ERROR: ERROR: cannot open lasreadertxt with file name 'C:\Users\jwiesehahn\Downloads\temp'

I guess this is because attribute href in the VPC does not point correctly to the referenced files?!

{
  "type": "FeatureCollection",
  "features": [
  {
    "type": "Feature",
    "stac_version": "1.0.0",
    "stac_extensions": [
      "https://stac-extensions.github.io/pointcloud/v1.0.0/schema.json",
      "https://stac-extensions.github.io/projection/v1.1.0/schema.json"
     ],
    "id": "MixedConifer",
    "geometry": {
      "coordinates": [
        [ [-111.204028964,34.457659149,0], [-111.203051179,34.457659149,0], [-111.203051179,34.458471515,0], [-111.204028964,34.458471515,0], [-111.204028964,34.457659149,0] ]
      ],
      "type": "Polygon"
    },
    "bbox": [-111.204028964, 34.457659149, 0, -111.203051179, 34.458471515, 0],
    "properties": {
      "datetime": "0-01-01T00:00:00Z",
      "pc:count": 37657,
      "pc:type": "lidar",
      "proj:bbox": [481260.000, 3812921.090, 481349.990, 3813010.990],
      "proj:wtk2": "PROJCRS[\"NAD83 / UTM zone 12N\",BASEGEOGCRS[\"NAD83\",DATUM[\"North American Datum 1983\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4269]],CONVERSION[\"UTM zone 12N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-111,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Engineering survey, topographic mapping.\"],AREA[\"North America - between 114°W and 108°W - onshore and offshore. Canada - Alberta; Northwest Territories; Nunavut; Saskatchewan. United States (USA) - Arizona; Colorado; Idaho; Montana; New Mexico; Utah; Wyoming.\"],BBOX[31.33,-114,84,-108]],ID[\"EPSG\",26912]]",
      "proj:epsg": 26912,
      "index:indexed": false
    },
    "links": [],
    "assets": {
      "data": {
      "href": "./",
      "roles": ["data"]
      }
    }
  }
  ]
}

did this change? The path looked different in https://github.com/r-lidar/lasR/issues/34#issue-2268660998

[Jean-Romain]

I cannot reproduce. With o = tempfile(fileext = ".vpc") I get a correct output. I guess it is a Windows specific issue. Yet this very specific uses cases it tested on windows in the unit tests here

{
  "type": "FeatureCollection",
  "features": [
  {
    "type": "Feature",
    "stac_version": "1.0.0",
    "stac_extensions": [
      "https://stac-extensions.github.io/pointcloud/v1.0.0/schema.json",
      "https://stac-extensions.github.io/projection/v1.1.0/schema.json"
     ],
    "id": "MixedConifer",
    "geometry": {
      "coordinates": [
        [ [34.457661675,-111.204027561,0], [34.458474039,-111.204027561,0], [34.458474039,-111.203049777,0], [34.457661675,-111.203049777,0], [34.457661675,-111.204027561,0] ]
      ],
      "type": "Polygon"
    },
    "bbox": [34.457661675, -111.204027561, 0, 34.458474039, -111.203049777, 0],
    "properties": {
      "datetime": "0-01-01T00:00:00Z",
      "pc:count": 37657,
      "pc:type": "lidar",
      "proj:bbox": [481260.000, 3812921.090, 481349.990, 3813010.990],
      "proj:wtk2": "PROJCRS[\"NAD83 / UTM zone 12N\",BASEGEOGCRS[\"NAD83\",DATUM[\"North American Datum 1983\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4269]],CONVERSION[\"UTM zone 12N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-111,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Engineering survey, topographic mapping.\"],AREA[\"North America - between 114°W and 108°W - onshore and offshore. Canada - Alberta; Northwest Territories; Nunavut; Saskatchewan. United States (USA) - Arizona; Colorado; Idaho; Montana; New Mexico; Utah; Wyoming.\"],BBOX[31.33,-114,84,-108]],ID[\"EPSG\",26912]]",
      "proj:epsg": 26912,
      "index:indexed": false
    },
    "links": [],
    "assets": {
      "data": {
      "href": "./../../home/jr/R/x86_64-pc-linux-gnu-library/4.4/lasR/extdata/MixedConifer.las",
      "roles": ["data"]
      }
    }
  }
  ]
}

[wiesehahn]

with o = tempfile(fileext = ".vpc") I also get

{
  "type": "FeatureCollection",
  "features": [
  {
    "type": "Feature",
    "stac_version": "1.0.0",
    "stac_extensions": [
      "https://stac-extensions.github.io/pointcloud/v1.0.0/schema.json",
      "https://stac-extensions.github.io/projection/v1.1.0/schema.json"
     ],
    "id": "MixedConifer",
    "geometry": {
      "coordinates": [
        [ [-111.204028964,34.457659149,0], [-111.203051179,34.457659149,0], [-111.203051179,34.458471515,0], [-111.204028964,34.458471515,0], [-111.204028964,34.457659149,0] ]
      ],
      "type": "Polygon"
    },
    "bbox": [-111.204028964, 34.457659149, 0, -111.203051179, 34.458471515, 0],
    "properties": {
      "datetime": "0-01-01T00:00:00Z",
      "pc:count": 37657,
      "pc:type": "lidar",
      "proj:bbox": [481260.000, 3812921.090, 481349.990, 3813010.990],
      "proj:wtk2": "PROJCRS[\"NAD83 / UTM zone 12N\",BASEGEOGCRS[\"NAD83\",DATUM[\"North American Datum 1983\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4269]],CONVERSION[\"UTM zone 12N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-111,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Engineering survey, topographic mapping.\"],AREA[\"North America - between 114°W and 108°W - onshore and offshore. Canada - Alberta; Northwest Territories; Nunavut; Saskatchewan. United States (USA) - Arizona; Colorado; Idaho; Montana; New Mexico; Utah; Wyoming.\"],BBOX[31.33,-114,84,-108]],ID[\"EPSG\",26912]]",
      "proj:epsg": 26912,
      "index:indexed": false
    },
    "links": [],
    "assets": {
      "data": {
      "href": "./",
      "roles": ["data"]
      }
    }
  }
  ]
}

locally the test creates a VPC with href as "href": "./testvpc/bcts_1.laz" ... in which case

pipeline = reader_las() + hulls()
            ans = exec(pipeline, on = ans)

runs fine as the VPC is in the parent folder of testvpc.

But runnng this for f <- system.file("extdata", "MixedConifer.las", package="lasR") creates a VPC with href as "href": "./" for me which gives the error.

---

Apart from href what is also strange is that lon and lat are changed under bbox and coordinates in your VPC compared to my output?!

[Jean-Romain]

Apart from href what is also strange is that lon and lat are changed under bbox and coordinates in your VPC compared to my output?!

That sucks! The transformation is handled by GDAL and I have no idea why the output are different. I will write more unit tests and run them with github action on all plateforms.

[Jean-Romain]

This is the file I just created on a windows machine. long lat are inverted but href is correct

{
  "type": "FeatureCollection",
  "features": [
  {
    "type": "Feature",
    "stac_version": "1.0.0",
    "stac_extensions": [
      "https://stac-extensions.github.io/pointcloud/v1.0.0/schema.json",
      "https://stac-extensions.github.io/projection/v1.1.0/schema.json"
     ],

    "id": "MixedConifer",
    "geometry": {
      "coordinates": [
        [ [-111.204028964,34.457659149,0], [-111.203051179,34.457659149,0], [-111.203051179,34.458471515,0], [-111.204028964,34.458471515,0], [-111.204028964,34.457659149,0] ]
      ],
      "type": "Polygon"
    },
    "bbox": [-111.204028964, 34.457659149, 0, -111.203051179, 34.458471515, 0],
    "properties": {
      "datetime": "0-01-01T00:00:00Z",
      "pc:count": 37657,
      "pc:type": "lidar",
      "proj:bbox": [481260.000, 3812921.090, 481349.990, 3813010.990],
      "proj:wtk2": "PROJCRS[\"NAD83 / UTM zone 12N\",BASEGEOGCRS[\"NAD83\",DATUM[\"North American Datum 1983\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4269]],CONVERSION[\"UTM zone 12N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-111,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Engineering survey, topographic mapping.\"],AREA[\"North America - between 114°W and 108°W - onshore and offshore. Canada - Alberta; Northwest Territories; Nunavut; Saskatchewan. United States (USA) - Arizona; Colorado; Idaho; Montana; New Mexico; Utah; Wyoming.\"],BBOX[31.33,-114,84,-108]],ID[\"EPSG\",26912]]",
      "proj:epsg": 26912,
      "index:indexed": false
    },
    "links": [],
    "assets": {
      "data": {
      "href": "./../../R/win-library/4.4/lasR/extdata/MixedConifer.las",
      "roles": ["data"]
      }
    }
  }
  ]
}

[wiesehahn]

I am a little confused, I thought lat/lon would be the standard for e.g. geojson, which in turn is mentioned for the VPC spec (https://github.com/PDAL/wrench/blob/main/vpc-spec.md#boundaries). But actually I think lon/lat is correct as this seems also to be the case in https://github.com/stac-extensions/pointcloud/blob/main/examples/example-autzen.json#L7 and QGIS reads it correctly like this.

So results under windows (here and in my case) seem to be correct.

About href I dont get yet why it is populated as "href": "./" in my case. Possibly something to do from which workspace it is created!?

It might make sense anyway to give the option to populate href either as relative or absolute path.

[Jean-Romain]

It might make sense anyway to give the option to populate href either as relative or absolute path.

Does it? Does Virtual Raster Tile supports absolute path? I have never created VRT raster with absolute path. PDAL creates VPC files with relative path. Does QGIS support absolute path? I'm not sure what the standard says about that.

I am a little confused, I thought lat/lon would be the standard for e.g. geojson, which in turn is mentioned for the VPC spec. But actually I think lon/lat is correct as this seems also to be the case in https://github.com/stac-extensions/pointcloud/blob/main/examples/example-autzen.json#L7 and QGIS reads it correctly like this.

Actually I copied the order you provided few day ago in the example file you created with QGIS. So it means that the order is incorrect on Linux at least. Maybe macOS.

About href I dont get yet why it is populated as "href": "./" in my case. Possibly something to do from which workspace it is created!?

This one requires investigation. Sadly I tested on Windows successfully.

[wiesehahn]

As VPC is a STAC catalog which is mainly used to reference data in cloud storage these kind of paths (not sure if S3 counts as absolute but I think so) are common (e.g. https://github.com/stac-extensions/pointcloud/blob/main/examples/example-autzen.json#L4).

Not sure about VRT, I thought it was possible but can not find it currently, but at least with gdaltindex (https://gdal.org/programs/gdaltindex.html#cmdoption-gdaltindex-write_absolute_path) and gti driver it is possible. I guess relative would be default though.

[Jean-Romain]

I pushed a commit with a new argument absolute_path for write_vpc

[Jean-Romain]

So the inversion of long lat is weird and critical. If I run write_vpc in a fresh R session I get the coordinates in the correct order. But if run a second time the function the coordinates get inverted. Can you confirm if you observe something similar?

I did not reproduce the href issue.

[wiesehahn]

So the inversion of long lat is weird and critical. If I run write_vpc in a fresh R session I get the coordinates in the correct order. But if run a second time the function the coordinates get inverted. Can you confirm if you observe something similar?

No, for me it is always in lon / lat order, didn't experience something else so far.

I did not reproduce the href issue.

I think the issue with href is related to the case when I try to write the VPC on a different drive than the data its pointing to. In my case I have a RStudio project with renv on some network drive, so in this case the data is also at e.g. "Y:/Username/project/renv/library/R-4.3/x86_64-w64-mingw32/lasR/extdata/MixedConifer.las".

When writing VPC to disk with

f <- system.file("extdata", "MixedConifer.las", package="lasR")
# "Y:/Username/project/renv/library/R-4.3/x86_64-w64-mingw32/lasR/extdata/MixedConifer.las"

o = here::here("data/test/test.vpc")
# "Y:/Username/project/data/test/test.vpc"

pipeline = reader_las() + write_vpc(o)
ans = exec(pipeline, on = f)

attribute href is "href": "./../../renv/library/R-4.3/x86_64-w64-mingw32/lasR/extdata/MixedConifer.las",

When writing VPC to disk with

f <- system.file("extdata", "MixedConifer.las", package="lasR")
# "Y:/Username/project/renv/library/R-4.3/x86_64-w64-mingw32/lasR/extdata/MixedConifer.las"

o = tempfile(fileext = ".vpc")
# "C:\\Users\\JWIESE~1\\AppData\\Local\\Temp\\RtmpquD1Mm\\file1a14208469a5.vpc"

pipeline = reader_las() + write_vpc(o)
ans = exec(pipeline, on = f)

attribute href is "href": "./"

[Jean-Romain]

I think the issue with href is related to the case when I try to write the VPC on a different drive than the data its pointing to.

Ha ok!! It makes sense. What is supposed to be the relative path from a drive to another on Windows? On Windows C: and D: do not have a common ancestor. A relative path between the two drives cannot exist. On Linux and Mac it should work since there is no concept of "drives" and everything have a common ancestor /

I tested this in C++ (it corresponds to internal code of lasR)

#include <iostream>
#include <filesystem>
#include <string>

namespace fs = std::filesystem;

int main() 
{
    fs::path file = "D:/directory1/file.txt";
    fs::path output_path = "C:/directory2/directory3";
    fs::path relative = fs::relative(file, output_path);
    std::string relative_path = relative.string();
    std::cout << "Relative path: " << relative_path << std::endl;
    return 0;
}

Surprisingly it returns ../../../D:/directory1/file.txt. But this is on Linux. I'm pretty sure it returns an empty string on Windows according to your tests. I don't think this output is a valid relative path on Windows. See also this SO question.

I added an absolute_path argument. It should cover your case. Also I added a test. If the relative path is empty it writes the absolute path. Let me know if it behaves properly on your machine.

[wiesehahn]

Yes, writing absolute path works automatic and with setting.