Missing GeoKey for horizontal datum when using custom LCC

[oregonparks]

I just ran into a problem after using las2las (version 120618) to a dataset from UTM z10 nad83 to a custom LCC nad83. The points were in the correct spot after the project operation but I noticed that the GeoKey tags in the las header didn't include a horizontal datum, only an ellipsoid. The result was that ArcGIS tools for Las Datasets and Mosaic Datasets did not correctly interpret the las file Spatial Reference and thus rendered the data in the wrong location.

Digging a little more it appears that the critical missing key is: GeogGeodeticDatumGeoKey, Key ID = 2050

My solution was to use the -remove_all_vlrs flag to get rid of the SR altogether. Not my preference but worked for this case.

[Martin Isenburg]

Hello,

can anyone confirm that this particular "GeogGeodeticDatumGeoKey" GeoTIFF key really is needed by ArcGIS? Below is the result produced by LAStools. Is the "GeogEllipsoidGeoKey" GeoTIFF key not the one that should be taking care of this? Below the example and attached a tiny sample LAZ file with only 50000 points.

las2las.exe -i 44116H7108A.laz ^

  -utm 11N -nad83 ^

  -target_lcc 1312335.958005249 0 feet 41.75 -120.5 43.0 45.5 ^

  -target_feet -target_elevation_feet ^

   -o 44116H7108Aout.laz 

using projection UTM 'UTM 11 northern hemisphere'

using ellipsoid '11 - GRS 1980 (NAD-83,GDA-94) (6.37814e+006 0.00669438)'

using LCC target projection 'false east/north: 400000/0 [m], origin lat/ meridian long: 41.75/-120.5, parallel 1st/2nd: 43/45.5'

D:\LAStools\bin>lasinfo 44116H7108Aout.laz

lasinfo for 44116H7108Aout.laz

reporting all LAS header entries:

  file signature:             'LASF'

  file source ID:             0

  global_encoding:            0

  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000

  version major.minor:        1.2

  system identifier:          'LAStools (c) by rapidlasso GmbH'

  generating software:        'las2las (version 150304)'

  file creation day/year:     0/0

  header size:                227

  offset to point data:       489

  number var. length records: 2

  point data format:          1

  point data record length:   28

  number of point records:    211471

  number of points by return: 211468 3 0 0 0

  scale factor x y z:         0.01 0.01 0.01

  offset x y z:               2200000 1200000 0

  min x y z:                  2256283.72 1199864.59 2337.60

  max x y z:                  2256826.97 1201093.87 2598.85

variable length header record 1 of 2:

  reserved             43707

  user ID              'LASF_Projection'

  record ID            34735

  length after header  104

  description          'by LAStools of rapidlasso GmbH'

    GeoKeyDirectoryTag version 1.1.0 number of keys 12

      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected

      key 3072 tiff_tag_location 0 count 1 value_offset 32767 - ProjectedCSTypeGeoKey: user-defined

      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 9002 - ProjLinearUnitsGeoKey: Linear_Foot

      key 3078 tiff_tag_location 34736 count 1 value_offset 0 - ProjStdParallel1GeoKey: 43

      key 3079 tiff_tag_location 34736 count 1 value_offset 1 - ProjStdParallel2GeoKey: 45.5

      key 3088 tiff_tag_location 34736 count 1 value_offset 2 - ProjCenterLongGeoKey: -120.5

      key 3081 tiff_tag_location 34736 count 1 value_offset 3 - ProjNatOriginLatGeoKey: 41.75

      key 3082 tiff_tag_location 34736 count 1 value_offset 4 - ProjFalseEastingGeoKey: 1312335.958

      key 3083 tiff_tag_location 34736 count 1 value_offset 5 - ProjFalseNorthingGeoKey: 0

      key 2056 tiff_tag_location 0 count 1 value_offset 7019 - GeogEllipsoidGeoKey: Ellipse_GRS_1980

      key 4099 tiff_tag_location 0 count 1 value_offset 9002 - VerticalUnitsGeoKey: Linear_Foot

variable length header record 2 of 2:

  reserved             43707

  user ID              'LASF_Projection'

  record ID            34736

  length after header  48

  description          'by LAStools of rapidlasso GmbH'

    GeoDoubleParamsTag (number of doubles 6)

      43 45.5 -120.5 41.75 1.31234e+006 0

the header is followed by 2 user-defined bytes

LASzip compression (version 2.2r0 c2 50000): POINT10 2 GPSTIME11 2

reporting minimum and maximum for all LAS point record entries ...

  X             5628372    5682697

  Y              -13541     109387

  Z              233760     259885

  intensity           1        222

  return_number       1          2

  number_of_returns   1          2

  edge_of_flight_line 0          0

  scan_direction_flag 0          1

  classification      1          8

  scan_angle_rank    -2         15

  user_data         115        132

  point_source_ID 15834      24836

  gps_time 421201.257315 422107.774912

overview over number of returns of given pulse: 211463 8 0 0 0 0 0

histogram of classification of points:

          139963  unclassified (1)

           58106  ground (2)

           13402  keypoint (8)

--
Download LAStools at
http://lastools.org
http://rapidlasso.com
Be social with LAStools at
http://facebook.com/LAStools
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe

[Paul Beaty]

Hello Martin,

I cannot open the file in ArgGIS for Desktop 10.3.

 

I can open in ERDAS IMAGINE 2015, but an error message “Unknown datum -1” pops up; with the session log stating the following:

07/03/15 10:27:09 pointcloudinfo -image C:/Temp/LiDAR/44116H7108Aout_sub.laz;

07/03/15 10:27:10 SessionMgr(5792): WARNING: #7931 from gk_CSInfoGCSReadUserDefined

07/03/15 10:27:10 SessionMgr(5792): WARNING: Unknown datum: -1

 

Paul Beaty

[Martin Isenburg]

Hi Paul,

Unlike Hexagon (Leica/Intergraph/ERDAS) who - along my many other companies (see http://laszip.org for a complete list) - support compressed LAZ in many of their products, the ESRI Empire - despite loud pleas from its user base - has so far resisted adding LASzip to its Arghhh ... (-:

After decompressing LAZ to LAS you should be able to import the file, no?

Does the "WARNING: Unknown datum: -1" message indicate that it is really the  "GeogGeodeticDatumGeoKey" GeoTIFF key that is missing ... ?

Martin

[Martin Isenburg]

Hello,

okay, I added the "GeogGeodeticDatumGeoKey" to the GeoTIFF tags. Does it work now? Small sample is attached again.

Regards,

Martin @rapidlasso

======================

D:\LAStools\bin>lasinfo 44116H7108Aout_sub2.laz

lasinfo for 44116H7108Aout_sub2.laz

reporting all LAS header entries:

  file signature:             'LASF'

  file source ID:             0

  global_encoding:            0

  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000

  version major.minor:        1.2

  system identifier:          'LAStools (c) by rapidlasso GmbH'

  generating software:        'las2las (version 150304)'

  file creation day/year:     0/0

  header size:                227

  offset to point data:       497

  number var. length records: 2

  point data format:          1

  point data record length:   28

  number of point records:    20000

  number of points by return: 20000 0 0 0 0

  scale factor x y z:         0.01 0.01 0.01

  offset x y z:               2200000 1200000 0

  min x y z:                  2256312.36 1199866.68 2485.73

  max x y z:                  2256559.39 1200432.72 2556.66

variable length header record 1 of 2:

  reserved             43707

  user ID              'LASF_Projection'

  record ID            34735

  length after header  112

  description          'by LAStools of rapidlasso GmbH'

    GeoKeyDirectoryTag version 1.1.0 number of keys 13

      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected

      key 3072 tiff_tag_location 0 count 1 value_offset 32767 - ProjectedCSTypeGeoKey: user-defined

      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 9002 - ProjLinearUnitsGeoKey: Linear_Foot

      key 3078 tiff_tag_location 34736 count 1 value_offset 0 - ProjStdParallel1GeoKey: 43

      key 3079 tiff_tag_location 34736 count 1 value_offset 1 - ProjStdParallel2GeoKey: 45.5

      key 3088 tiff_tag_location 34736 count 1 value_offset 2 - ProjCenterLongGeoKey: -120.5

      key 3081 tiff_tag_location 34736 count 1 value_offset 3 - ProjNatOriginLatGeoKey: 41.75

      key 3082 tiff_tag_location 34736 count 1 value_offset 4 - ProjFalseEastingGeoKey: 1312335.958

      key 3083 tiff_tag_location 34736 count 1 value_offset 5 - ProjFalseNorthingGeoKey: 0

      key 2050 tiff_tag_location 0 count 1 value_offset 6019 - GeogGeodeticDatumGeoKey: DatumE_GRS1980

      key 2056 tiff_tag_location 0 count 1 value_offset 7019 - GeogEllipsoidGeoKey: Ellipse_GRS_1980

      key 4099 tiff_tag_location 0 count 1 value_offset 9002 - VerticalUnitsGeoKey: Linear_Foot

variable length header record 2 of 2:

  reserved             43707

  user ID              'LASF_Projection'

  record ID            34736

  length after header  48

  description          'by LAStools of rapidlasso GmbH'

    GeoDoubleParamsTag (number of doubles 6)

      43 45.5 -120.5 41.75 1.31234e+006 0

the header is followed by 2 user-defined bytes

LASzip compression (version 2.2r0 c2 50000): POINT10 2 GPSTIME11 2

reporting minimum and maximum for all LAS point record entries ...

  X             5631236    5655939

  Y              -13332      43272

  Z              248573     255666

  intensity          59        209

  return_number       1          1

  number_of_returns   1          1

  edge_of_flight_line 0          0

  scan_direction_flag 0          1

  classification      1          8

  scan_angle_rank     4         15

  user_data         118        130

  point_source_ID 15834      24833

  gps_time 421201.257315 421204.516034

overview over number of returns of given pulse: 20000 0 0 0 0 0 0

histogram of classification of points:

           13464  unclassified (1)

            5746  ground (2)

             790  keypoint (8)

[oregonparks]

Almost there. In your second example the GeoKey is recognized by ArcGIS now. See attached 'YesDatumKey'

However, the particular datum choice is not correct. Your second sample uses an Ellipsoid Only datum (DatumE_XXX) of GRS80 when it should be 'Datum_North_American_Datum_1983' since the original data was 'UTM 11n -nad83' which would reference a Geodetic Datum of NAD83 as opposed to an Ellipsoid.    http://www.remotesensing.org/geotiff/spec/geotiff6.html#6.3.2.2

See the attached GeoKey_Difference.jpg for the results.

Regards,

Brady

[Martin Isenburg]

Hello,

hope it's correct now. Where does this difference come from? What makes a the specification of the horizontal datum NAD83 difference from the ellipsoid that this horizontal datum's LCC projection is defined on?

Martin

>>las2las.exe -i 44116H7108A.laz ^

                       -utm 11N -nad83 ^

                       -target_lcc 1312335.958005249 0 feet 41.75 -120.5 43.0 45.5 ^

                       -target_feet -target_elevation_feet ^

                       -subseq 0 30000 ^

                       -o 44116H7108Aout_sub3.laz 

using projection UTM 'UTM 11 northern hemisphere'

using ellipsoid '11 - GRS 1980 (NAD-83,GDA-94) (6.37814e+006 0.00669438)'

using LCC target projection 'false east/north: 400000/0 [m], origin lat/ meridian long: 41.75/-120.5, parallel 1st/2nd: 43/45.5'

>> lasinfo 44116H7108Aout_sub3.laz

lasinfo for 44116H7108Aout_sub3.laz

reporting all LAS header entries:

  file signature:             'LASF'

  file source ID:             0

  global_encoding:            0

  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000

  version major.minor:        1.2

  system identifier:          'LAStools (c) by rapidlasso GmbH'

  generating software:        'las2las (version 150304)'

  file creation day/year:     0/0

  header size:                227

  offset to point data:       497

  number var. length records: 2

  point data format:          1

  point data record length:   28

  number of point records:    30000

  number of points by return: 30000 0 0 0 0

  scale factor x y z:         0.01 0.01 0.01

  offset x y z:               2200000 1200000 0

  min x y z:                  2256307.20 1199866.68 2473.36

  max x y z:                  2256607.47 1200551.72 2574.18

variable length header record 1 of 2:

  reserved             43707

  user ID              'LASF_Projection'

  record ID            34735

  length after header  112

  description          'by LAStools of rapidlasso GmbH'

    GeoKeyDirectoryTag version 1.1.0 number of keys 13

      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected

      key 3072 tiff_tag_location 0 count 1 value_offset 32767 - ProjectedCSTypeGeoKey: user-defined

      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 9002 - ProjLinearUnitsGeoKey: Linear_Foot

      key 3078 tiff_tag_location 34736 count 1 value_offset 0 - ProjStdParallel1GeoKey: 43

      key 3079 tiff_tag_location 34736 count 1 value_offset 1 - ProjStdParallel2GeoKey: 45.5

      key 3088 tiff_tag_location 34736 count 1 value_offset 2 - ProjCenterLongGeoKey: -120.5

      key 3081 tiff_tag_location 34736 count 1 value_offset 3 - ProjNatOriginLatGeoKey: 41.75

      key 3082 tiff_tag_location 34736 count 1 value_offset 4 - ProjFalseEastingGeoKey: 1312335.958

      key 3083 tiff_tag_location 34736 count 1 value_offset 5 - ProjFalseNorthingGeoKey: 0

      key 2050 tiff_tag_location 0 count 1 value_offset 6269 - GeogGeodeticDatumGeoKey: Datum_North_American_Datum_1983

      key 2056 tiff_tag_location 0 count 1 value_offset 7019 - GeogEllipsoidGeoKey: Ellipse_GRS_1980

      key 4099 tiff_tag_location 0 count 1 value_offset 9002 - VerticalUnitsGeoKey: Linear_Foot

variable length header record 2 of 2:

  reserved             43707

  user ID              'LASF_Projection'

  record ID            34736

  length after header  48

  description          'by LAStools of rapidlasso GmbH'

    GeoDoubleParamsTag (number of doubles 6)

      43 45.5 -120.5 41.75 1.31234e+006 0

the header is followed by 2 user-defined bytes

LASzip compression (version 2.2r0 c2 50000): POINT10 2 GPSTIME11 2

reporting minimum and maximum for all LAS point record entries ...

  X             5630720    5660747

  Y              -13332      55172

  Z              247336     257418

  intensity          59        209

  return_number       1          1

  number_of_returns   1          1

  edge_of_flight_line 0          0

  scan_direction_flag 0          1

  classification      1          8

  scan_angle_rank     4         15

  user_data         118        130

  point_source_ID 15834      24833

  gps_time 421201.257315 421205.196815

overview over number of returns of given pulse: 30000 0 0 0 0 0 0

histogram of classification of points:

           21266  unclassified (1)

            7544  ground (2)

            1190  keypoint (8)

[Evon Silvia]

If I understand the problem correction, the difference is that the ellipsoid is simply a mathematical shape defined by two radii or one radius and a flattening ratio, while the datum defines the position of that ellipsoid relative to the surface and/or center of the Earth. Although it is true that every datum must have one reference ellipsoid, it is not necessarily true that every ellipsoid only has one datum - i.e., there's not necessarily a 1-to-1 relationship between ellipsoids and horizontal datums.

Take the following graphic with an ellipsoid as defined by its major and minor radii and an exaggeration of Earth's irregular surface. Below that I show three possible configurations of the exact same ellipsoid relative to the earth - a bad fit, a biased fit, and a best fit. There are many more, but these three make the distinction I'm attempting.

Although the same ellipsoid is used every time, the way that it fits the center and surface of the Earth varies. This is implicitly defined by the datum definition. 

1. The bad fit is mostly shown as an example, but I'm sure that someone can come up with a scenario where it makes sense - maybe an extremely localized datum?
2. The biased fit represents an attempt to match the ellipsoid to a particular region; this is the approach taken for NAD83, which tries to fix the datum to the North American tectonic plate.
3. The best fit, by contrast, attempts to match the ellipsoid to the center (or surface) of the earth as closely as possible; this is the approach taken for WGS84.
   

Therefore, defining the ellipsoid does not provide enough information to uniquely define your file's coordinate system. Thus the need to also define the horizontal datum.

Hope that helps and I didn't screw that up too much this early on a Monday morning. :) I'll let the other geodesists on this list correct me as needed.

Evon

--

Evon Silvia PLS

Geomatics Specialist

Quantum Spatial

517 SW 2nd Street, Suite 400, Corvallis, OR 97333

P: (541) 752-1204

E: esi...@quantumspatial.com

[oregonparks]

Thanks Evon for the cool graphic and explanation of how the ellipsoid and the datum work together!

Martin's 44116H7108Aout_sub3.laz sample seemed to be correct at first but then I looked at the coordinate system definition of the Las Dataset in ArcMap and noticed that the semiminor axis and inverse flattening were incorrectly defined for GRS_80. See attached.

Perhaps Martin's question about this particular GeoKey is a good one. On review, there is apparently another key that may be relevant: GeographicTypeGeoKey, Key ID = 2048. The name of this key doesn't seem quite correct but the possible values include a GCS_NAD83.

I contacted ESRI's experts in both Coordinate Systems and LiDAR and sadly neither one knew how the GeoTIFF GeoKeys are mapped to an ArcGIS coordinate system definition.

Martin, would it be possible to add GeographicTypeGeoKey = GCS_NAD83  instead of the GeogGeodeticDatumGeoKey, Key ID = 2050 and see if that is correct?