search

r-lidar / rlas

R package to read and write las and laz files used to store LiDAR data
https://cran.r-project.org/package=rlas
GNU General Public License v3.0
34 stars 14 forks source link

CRS of las/laz file #1

Closed edzer closed 7 years ago

edzer commented 7 years ago

Does rlas give access to the coordinate reference system of the las/laz file? I read that version 1.4 should support this. Am asking because it would be nice to set up exchange with simple features in R, and through that to leafet/mapview or 3D visualisations, the kind of things @mdsumner and @tim-salabim do.

See also https://github.com/tim-salabim/rspatial/issues/6

Jean-Romain commented 7 years ago

You're right version 1.4 should support that information. But:

  1. I have never encounter such file so I cannot make tests
  2. I'm not comfortable with spatial data, spatial projections, CRS and so on. It's out of my knowledge.

According to point 2 I don't really understand what do you want and I don't really understand what are the purposes of the your package or the one written by mdsummer. Are you expected a feature to switch the CRS from UTM to NAD for example?

mdsumner commented 7 years ago

@Jean-Romain it's really for being able to map the coordinates with other data, doing anything that requires coincidence It's trivial to transform the values using functions like rgdal::project and now sf::st_transform, but the source projection (family, central longitude, latitude, datum, ellipsoid, false offsets, scalings and others) is critical to doing that.

It should be relatively simple if it's in the file, a PROJ.4 string, or an EPSG code, or a table of parameters. Most LiDAR files I've seen come with an assumption about using the "local system", which is usually UTM in some form.

I believe older forms used GeoTIFF tags for this information. and in my lastools (that is pretty out of date) I see the output below, but I don't get the GeoTIFF parts in the output of rlas::readlasheader.

It's probably different for 1.4 files, but perhaps it's not too difficult to extract with the libs?

 ../bin/lasinfo test.laz
reporting all LAS header entries:
  file signature:             'LASF'
 ...
...
  length after header  192
  description          'GeoTiff Projection Keys'
    GeoKeyDirectoryTag version 1.1.0 number of keys 23
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 2048 tiff_tag_location 0 count 1 value_offset 4269 - GeographicTypeGeoKey: GCS_NAD83
      key 2049 tiff_tag_location 34737 count 24 value_offset 54 - GeogCitationGeoKey: GCS_North_American_1983
      key 2050 tiff_tag_location 0 count 1 value_offset 6269 - GeogGeodeticDatumGeoKey: Datum_North_American_Datum_1983
      key 2051 tiff_tag_location 0 count 1 value_offset 8901 - GeogPrimeMeridianGeoKey: PM_Greenwich
      key 2054 tiff_tag_location 0 count 1 value_offset 9102 - GeogAngularUnitsGeoKey: Angular_Degree
      key 2055 tiff_tag_location 34736 count 1 value_offset 11 - GeogAngularUnitSizeGeoKey: 0.01745329252
      key 2056 tiff_tag_location 0 count 1 value_offset 7019 - GeogEllipsoidGeoKey: Ellipse_GRS_1980
      key 2057 tiff_tag_location 34736 count 1 value_offset 8 - GeogSemiMajorAxisGeoKey: 6378137
      key 2059 tiff_tag_location 34736 count 1 value_offset 9 - GeogInvFlatteningGeoKey: 298.2572221
      key 2061 tiff_tag_location 34736 count 1 value_offset 10 - GeogPrimeMeridianLongGeoKey: 0
      key 3072 tiff_tag_location 0 count 1 value_offset 32128 - ProjectedCSTypeGeoKey: PCS_NAD83_Pennsylvania_N
      key 3073 tiff_tag_location 34737 count 54 value_offset 0 - PCSCitationGeoKey: NAD_1983_StatePlane_Pennsylvania_North_FIPS_3701_Feet
      key 3075 tiff_tag_location 0 count 1 value_offset 8 - ProjCoordTransGeoKey: CT_LambertConfConic_2SP
      key 3076 tiff_tag_location 0 count 1 value_offset 9003 - ProjLinearUnitsGeoKey: Linear_Foot_US_Survey
      key 3077 tiff_tag_location 34736 count 1 value_offset 7 - ProjLinearUnitSizeGeoKey: 0.3048006096
      key 3078 tiff_tag_location 34736 count 1 value_offset 3 - ProjStdParallel1GeoKey: 40.88333333
      key 3079 tiff_tag_location 34736 count 1 value_offset 4 - ProjStdParallel2GeoKey: 41.95
      key 3081 tiff_tag_location 34736 count 1 value_offset 6 - ProjNatOriginLatGeoKey: 40.16666667
      key 3082 tiff_tag_location 34736 count 1 value_offset 0 - ProjFalseEastingGeoKey: 1968500
      key 3083 tiff_tag_location 34736 count 1 value_offset 1 - ProjFalseNorthingGeoKey: 0
      key 3088 tiff_tag_location 34736 count 1 value_offset 2 - ProjCenterLongGeoKey: -77.75
      key 3093 tiff_tag_location 34736 count 1 value_offset 5 - ProjScaleAtCenterGeoKey: 0
variable length header record 2 of 5:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34736
  length after header  96
  description          'GeoTiff double parameters'
    GeoDoubleParamsTag (number of doubles 12)
      1.9685e+006 0 -77.75 40.8833 41.95 0 40.1667 0.304801 6.37814e+006 298.257 0 0.0174533
variable length header record 3 of 5:
edzer commented 7 years ago

Yes, the CRS allows you to understand where exactly on Earth the data were collected, and then to combine this with any other information for which you know.

The standard mentions it is (should) be encoded as WKT, which is essentially a text description that is understood by other software.

Jean-Romain commented 7 years ago

Ok, understood and I agree with you.

LAS specifications 1.0 to 1.3 do not have standard field to store the CRS. It is usually stored in a kind of extra data if I well understood the LAS specifications. But most of the time (in my own experience) there is no information relative to CRS. Anyway this information is an array of strings and it should not be difficult to retrieve these strings in the header if they exist. rlas rely on LASlib so this feature is already implemented. I simply forgot to wrap it into R. I will study the LASlib code to find where is the function to extract such data.

Jean-Romain commented 7 years ago

I figured out how to read these data. More difficult than expected (maybe this is why I forgot to wrap this feature :smile: ). I will need your help to understand what the files actually contain.

The files do not contain as much data as you copied in your post. Actually there is a post process on the raw data to get this output. It seems to contain mostly some code/tag/key that refer to some CRS and other spatial data. Here what I can easily extract from one of my files:

variable length header record 1 of 5
  reserved             43707
  user ID              LASF_Projection
  record ID            34735
  length after header  176
  description          GeoTiff Projection Keys
  data                 
variable length header record 2 of 5
  reserved             43707
  user ID              LASF_Projection
  record ID            34736
  length after header  80
  description          GeoTiff double parameters
  data                 
variable length header record 3 of 5
  reserved             43707
  user ID              LASF_Projection
  record ID            34737
  length after header  96
  description          GeoTiff ASCII parameters
  data                 NAD_1983_CSRS98_UTM_Zone_17N|GCS_North_American_1983_CSRS98|datum: D_North_American_1983_CSRS98|p
variable length header record 4 of 5
  reserved             43707
  user ID              NIIRS10
  record ID            4
  length after header  10
  description          NIIRS10 Timestamp
  data                 
variable length header record 5 of 5
  reserved             43707
  user ID              NIIRS10
  record ID            1
  length after header  26
  description          NIIRS10 Tile Index
  data

To go deeper in the content of the header I have to write many many conditional statements on the tags/keys to process the different flavor of data. That's what is done in lasinfo:

if (lasheader->vlrs[i].record_id == 34735) // GeoKeyDirectoryTag
{
  [...]
  switch (lasreader->header.vlr_geo_key_entries[j].key_id)
  {
    case 1024: // GTModelTypeGeoKey 
    switch (lasreader->header.vlr_geo_key_entries[j].value_offset)
    {
      case 1: // ModelTypeProjected   
        fprintf(file_out, "GTModelTypeGeoKey: ModelTypeProjected\012");
        break;
       case 2:
         fprintf(file_out, "GTModelTypeGeoKey: ModelTypeGeographic\012");
         break;
       [...]

So I think this is the deepest inspection I can do. To go further I have to copy/paste the whole set of conditional statement from lasinfo because I'm not able to write it myself.

I will work on that. But it will take me a while. The first level of conditional statement is the following one

if ((strcmp(lasheader->vlrs[i].user_id, "LASF_Projection") == 0) && (lasheader->vlrs[i].data != 0))
{
  if (lasheader->vlrs[i].record_id == 34735) // GeoKeyDirectoryTag
  {}
  else if (lasheader->vlrs[i].record_id == 34736) // GeoDoubleParamsTag
  {}
  else if (lasheader->vlrs[i].record_id == 34737) // GeoAsciiParamsTag
  {}
  else if (lasheader->vlrs[i].record_id == 2111) // WKT OGC MATH TRANSFORM
  {}
  else if (lasheader->vlrs[i].record_id == 2112) // WKT OGC COORDINATE SYSTEM
  {}
}
else if ((strcmp(lasheader->vlrs[i].user_id, "LASF_Spec") == 0) && (lasheader->vlrs[i].data != 0))
{}

What is the meaning of these tags GeoKeyDirectoryTag, GeoDoubleParamsTag, GeoAsciiParamsTag, WKT OGC MATH TRANSFORM and WKT OGC COORDINATE SYSTEM?

edzer commented 7 years ago

I don't think this will be very useful, since these GeoTIFF tags are not very interoperable (= understood by other software), if I understand it correctly. On page 4, the LAS 1.4 standard mentions that WKT (which is interoperable) will replace GeoTIFF tags for the CRS for point types 6-10. Maybe we should wait for some files in that format before we make a serious effort.

Jean-Romain commented 7 years ago

I trust you. I don't have any expertise on these questions. The only thing I know is that I'm likely to be able to read extra data and add them in the output of the readlasheader. It's a boring job so if you think it's not useful I'm not going to do it.

If you encounter a LAS 1.4 file please share it I will try to support the format.

edzer commented 7 years ago

I found one in the pdal.org here. When running there docker image,

docker run -v `pwd`:/data pdal/pdal:1.4 pdal info --all /data/autzen.laz -p 0

gives a lot of stuff, including the SRS:

  {
    "comp_spatialreference": "PROJCS[\"NAD_1983_HARN_Lambert_Conformal_Conic\",GEOGCS[\"GCS_North_American_1983_HARN\",DATUM[\"NAD83_High_Accuracy_Reference_Network\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6152\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",43],PARAMETER[\"standard_parallel_2\",45.5],PARAMETER[\"latitude_of_origin\",41.75],PARAMETER[\"central_meridian\",-120.5],PARAMETER[\"false_easting\",1312335.958005249],PARAMETER[\"false_northing\",0],UNIT[\"foot\",0.3048,AUTHORITY[\"EPSG\",\"9002\"]]]",

but it's pretty messy - the CRS is repeated many times, later on.

mdsumner commented 7 years ago

Can we extract the header as text in raw form, perhaps as an internal R function? Even un-parsed it would be helpful and then it's an R-level job to parse and interpret it.

Jean-Romain commented 7 years ago

Yes, certainly. I going to work on that today.

mdsumner commented 7 years ago

Great, thanks!

Jean-Romain commented 7 years ago

@edzer the file you gave me in attachment is a LAS 1.2 with GeoTiff tags as the one I parsed in my previous message. The information you want is contained in the GeoTiff ASCII parameters tag (see my previous post). So I'm already able to extract it. I will be able to release a new version to read such information today (but not to write it).

edzer commented 7 years ago

OK, thanks, and apologies for the confusion! I must have mixed up the pdal version with the las version...

Jean-Romain commented 7 years ago

I updated the devel branch. Install the new version from this branch with

devtools::install_github("Jean-Romain/rlas", ref = "devel" )

It's surely not bug free because I didn't test it on many files (4 files from 4 different projects). I also added a print function for LASheader objects. For the file you gave me the output is given below. You can access to the GeoAsciiTag with:

h = readlasheader(...)
h$`Variable Length Records`[[x]][[6]]

with x being the position of this information into the list. It might be x = 3 but I suppose it can change depending on the data stored in the file.

There is currently no support for the function writelas. You will be able to read the metadata but not to write it.

> print(h)

File signature:           LASF 
File source ID:           0 
Global encoding:          0 
Project ID - GUID:        0 
Version:                  01 . 02 
System identifier:         
Generating software:      TerraScan 
File creation d/y:        198 / 2013 
header size:              227 
Offset to point data:     1390 
Num. var. length record:  4 
Point data format:        03 
Point data record length: 34 
Num. of point records:    10653336 
Num. of points by return: 9100899 1238442 283980 30015 0 
Scale factor X Y Z:       0.01 0.01 0.01 
Offset X Y Z:             635577.8 848882.2 406.14 
min X Y Z:                635577.8 848882.2 406.14 
max X Y Z:                639003.7 853537.7 615.26 
Variable length records: 
   Variable length records 1/4
       Reserved:             43707 
       User ID:              LASF_Projection 
       record ID:            34735 
       Length after heaader: 184 
       Description:          GeoTIFF GeoKeyDirectoryTag 
       Tags:
          Key 1024 tiff_tag_location 0 count 1 value offset 1 
          Key 1025 tiff_tag_location 0 count 1 value offset 1 
          Key 1026 tiff_tag_location 34737 count 38 value offset 0 
          Key 2048 tiff_tag_location 0 count 1 value offset 32767 
          Key 2049 tiff_tag_location 34737 count 60 value offset 38 
          Key 2050 tiff_tag_location 0 count 1 value offset 6152 
          Key 2054 tiff_tag_location 0 count 1 value offset 9102 
          Key 2057 tiff_tag_location 34736 count 1 value offset 7 
          Key 2059 tiff_tag_location 34736 count 1 value offset 6 
          Key 2061 tiff_tag_location 34736 count 1 value offset 8 
          Key 3059 tiff_tag_location 0 count 1 value offset 1 
          Key 3072 tiff_tag_location 0 count 1 value offset 32767 
          Key 3074 tiff_tag_location 0 count 1 value offset 32767 
          Key 3075 tiff_tag_location 0 count 1 value offset 8 
          Key 3076 tiff_tag_location 0 count 1 value offset 9002 
          Key 3078 tiff_tag_location 34736 count 1 value offset 2 
          Key 3079 tiff_tag_location 34736 count 1 value offset 3 
          Key 3084 tiff_tag_location 34736 count 1 value offset 1 
          Key 3085 tiff_tag_location 34736 count 1 value offset 0 
          Key 3086 tiff_tag_location 34736 count 1 value offset 4 
          Key 3087 tiff_tag_location 34736 count 1 value offset 5 
   Variable length records 2/4
       Reserved:             43707 
       User ID:              LASF_Projection 
       record ID:            34736 
       Length after heaader: 72 
       Description:          GeoTIFF GeoDoubleParamsTag 
       data:                 41.75 -120.5 43 45.5 1312336 0 298.2572 6378137 0 
   Variable length records 3/4
       Reserved:             43707 
       User ID:              LASF_Projection 
       record ID:            34737 
       Length after heaader: 99 
       Description:          GeoTIFF GeoAsciiParamsTag 
       data:                 NAD_1983_HARN_Lambert_Conformal_Conic|GCS Name = GCS_North_American_1983_HARN|Primem = Greenwich|| 
   Variable length records 4/4
       Reserved:             43707 
       User ID:              liblas 
       record ID:            2112 
       Length after heaader: 592 
       Description:          OGR variant of OpenGIS WKT SRS
edzer commented 7 years ago

Thanks, I tried it for autzen.laz, but I see nothing in there (or in the output you give above) that can be parsed as a valid well-known-text representation of a CRS.

Jean-Romain commented 7 years ago

I can make nothing more for you in that case. This is the whole information contained in the file. It seems to me that the output given by pdal can be parsed from these data. But if the pdal output is not enough for you, it's a dead end issue.

Jean-Romain commented 7 years ago

Maybe that link can help you: http://gis.stackexchange.com/questions/173111/converting-geotiff-projection-definition-to-proj4

The whole list is now labelled:

h$`Variable Length Records`[[x]][[6]]

becomes

h$`Variable Length Records`$GeoAsciiParamsTag$tags