NAD83 (HARN) to WGS84

69 views
Skip to first unread message

Martin Isenburg

unread,
Feb 19, 2019, 2:23:39 AM2/19/19
to LAStools - efficient tools for LiDAR processing
Hello,

a user recently went to download LiDAR from the Illinois Height Modernization Program (aka ILHMP) and that acronym ILHMP pretty accurately transcribes the sounds that I made when I noticed that they offer uncompressed LAS of zipped folders of LAS for download... (-;

<ILHMP>

Because LASzip compression works really really well on these files:

LAS: 103,919,811 POINTS_112.las
LAZ:   14,611,758 POINTS_112.laz

LAS: 380,828,243 16SCJ1172.las
LAZ:   37,664,932 16SCJ1172.laz (*)

(*) after removing "fluff" reported by lasinfo by adding '-rescale 0.01 0.01 0.01' to the laszip command line

<ILHMP>

But at least the state of Illinois offers free and open LiDAR available via a nice portal (unlike most states in my native Germany who keep their LiDAR data locked up).

I've downloaded the two tiles shown in the two screenshots via the interactive portal. It states that the horizontal coordinate reference system is state plane "Illinois East" and "Illinois West" with NAD 1983 HARN as the horizontal datum and US survey feet as horizontal units. What calculation do I need to do to convert them into long/lat in the WGS84 datum (EPSG 4326). There were a number of previous discussions on the issue:


I'd like to go through those calculations myself. I understand there is no simple seven parameter Helmert transform that I could perform on the coordinates after transforming them to an ECEF representation such as 

las2las .... -transform_helmert -199.87,74.79,246.62 ...
las2las .... -transform_helmert 598.1,73.7,418.2,0.202,0.045,-2.455,6.7 ...

but instead there are correction grids per area that need to be used. Has someone seen where this is described step by step?

Regards,

Martin @rapidlasso

C:\software\LAStools\bin>lasinfo -i C:\Users\martin\Downloads\POINTS_112.las
lasinfo (190219) report for 'C:\Users\martin\Downloads\POINTS_112.las'
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   6C002954-AFBD-4F7E-5B9B-D15F73CD567F
  version major.minor:        1.2
  system identifier:          'NIIRS10'
  generating software:        'LIDAR1 tiled'
  file creation day/year:     358/2009
  header size:                227
  offset to point data:       2403
  number var. length records: 5
  point data format:          1
  point data record length:   28
  number of point records:    3711336
  number of points by return: 3193824 426909 86637 3966 0
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 0 0
  min x y z:                  1082336.00 1912295.00 619.04
  max x y z:                  1087336.00 1917294.99 942.26
variable length header record 1 of 5:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  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 32767 - GeographicTypeGeoKey: user-defined
      key 2049 tiff_tag_location 34737 count 63 value_offset 78 - GeogCitationGeoKey: GCS_North_American_1983_HARN|datum: D_North_American_1983_HARN
      key 2050 tiff_tag_location 0 count 1 value_offset 32767 - GeogGeodeticDatumGeoKey: user-defined
      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 9 - 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 6 - GeogSemiMajorAxisGeoKey: 6378137
      key 2059 tiff_tag_location 34736 count 1 value_offset 7 - GeogInvFlatteningGeoKey: 298.2572221
      key 2061 tiff_tag_location 34736 count 1 value_offset 8 - GeogPrimeMeridianLongGeoKey: 0
      key 3072 tiff_tag_location 0 count 1 value_offset 32767 - ProjectedCSTypeGeoKey: user-defined
      key 3073 tiff_tag_location 34737 count 54 value_offset 0 - PCSCitationGeoKey: NAD_1983_HARN_StatePlane_Illinois_East_FIPS_1201_Feet
      key 3075 tiff_tag_location 0 count 1 value_offset 1 - ProjCoordTransGeoKey: CT_TransverseMercator
      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 5 - ProjLinearUnitSizeGeoKey: 0.3048006096
      key 3081 tiff_tag_location 34736 count 1 value_offset 4 - ProjNatOriginLatGeoKey: 36.66666667
      key 3082 tiff_tag_location 34736 count 1 value_offset 0 - ProjFalseEastingGeoKey: 984250
      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: -88.33333333
      key 3092 tiff_tag_location 34736 count 1 value_offset 3 - ProjScaleAtNatOriginGeoKey: 0.999975
      key 4097 tiff_tag_location 34737 count 24 value_offset 54 - VerticalCitationGeoKey: NAVD88 - Geoid03 (Feet)
      key 4099 tiff_tag_location 0 count 1 value_offset 9003 - VerticalUnitsGeoKey: Linear_Foot_US_Survey
variable length header record 2 of 5:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34736
  length after header  80
  description          'GeoTiff double parameters'
    GeoDoubleParamsTag (number of doubles 10)
      984250 0 -88.3333 0.999975 36.6667 0.304801 6.37814e+006 298.257 0 0.0174533
variable length header record 3 of 5:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34737
  length after header  142
  description          'GeoTiff ASCII parameters'
    GeoAsciiParamsTag (number of characters 142)
      NAD_1983_HARN_StatePlane_Illinois_East_FIPS_1201_Feet|NAVD88 - Geoid03 (Feet)|GCS_North_American_1983_HARN|datum: D_North_American_1983_HARN|
variable length header record 4 of 5:
  reserved             43707
  user ID              'NIIRS10'
  record ID            4
  length after header  10
  description          'NIIRS10 Timestamp'
variable length header record 5 of 5:
  reserved             43707
  user ID              'NIIRS10'
  record ID            1
  length after header  26
  description          'NIIRS10 Tile Index'
the header is followed by 1456 user-defined bytes
reporting minimum and maximum for all LAS point record entries ...
  X           108233600  108733600
  Y           191229500  191729499
  Z               61904      94226
  intensity           1       5100
  return_number       1          4
  number_of_returns   1          4
  edge_of_flight_line 0          0
  scan_direction_flag 0          0
  classification      1         12
  scan_angle_rank     0          0
  user_data           0          0
  point_source_ID     8         11
  gps_time 377260.733428 379655.999837
number of first returns:        3193824
number of intermediate returns: 90679
number of last returns:         3193411
number of single returns:       2766578
overview over number of returns of given pulse: 2766578 680709 248172 15877 0 0 0
histogram of classification of points:
          924675  unclassified (1)
         1594586  ground (2)
             370  noise (7)
             287  water (9)
         1191418  overlap (12)
 +-> flagged as withheld:  1191418


C:\software\LAStools\bin>lasinfo -i C:\Users\martin\Downloads\16SCJ1172.las
lasinfo (190219) report for 'C:\Users\martin\Downloads\16SCJ1172.las'
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            1
  project ID GUID data 1-4:   7B0D6808-A5DC-4181-1081-5D81B4D8856B
  version major.minor:        1.2
  system identifier:          'ALS70'
  generating software:        'GeoCue LAS Updater'
  file creation day/year:     171/2016
  header size:                227
  offset to point data:       971
  number var. length records: 3
  point data format:          1
  point data record length:   28
  number of point records:    13600974
  number of points by return: 12831281 616941 142082 10458 212
  scale factor x y z:         0.001 0.001 0.001
  offset x y z:               2000000 950000 0
  min x y z:                  2568284.000 1023818.000 596.760
  max x y z:                  2573564.000 1029098.000 719.240
variable length header record 1 of 3:
  reserved             0
  user ID              'LASF_Projection'
  record ID            34735
  length after header  208
  description          'GeoTiff Projection Keys'
    GeoKeyDirectoryTag version 1.1.0 number of keys 25
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 1026 tiff_tag_location 34737 count 65 value_offset 86 - GTCitationGeoKey: PCS Name = NAD_1983_HARN_StatePlane_Illinois_West_FIPS_1202_Feet
      key 2048 tiff_tag_location 0 count 1 value_offset 4152 - GeographicTypeGeoKey: look-up for 4152 not implemented
      key 2049 tiff_tag_location 34737 count 113 value_offset 176 - GeogCitationGeoKey: GCS Name = GCS_North_American_1983_HARN|Datum = North_American_1983_HARN
|Ellipsoid = GRS_1980|Primem = Greenwich
      key 2050 tiff_tag_location 0 count 1 value_offset 6152 - GeogGeodeticDatumGeoKey: look-up for 6152 not implemented
      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 9 - 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 6 - GeogSemiMajorAxisGeoKey: 6378137
      key 2059 tiff_tag_location 34736 count 1 value_offset 7 - GeogInvFlatteningGeoKey: 298.2572221
      key 2061 tiff_tag_location 34736 count 1 value_offset 8 - GeogPrimeMeridianLongGeoKey: 0
      key 3059 tiff_tag_location 0 count 1 value_offset 1 - key ID 3059 not implemented
      key 3072 tiff_tag_location 0 count 1 value_offset 3444 - ProjectedCSTypeGeoKey: NAD83(HARN) / Illinois West (ftUS)
      key 3073 tiff_tag_location 34737 count 86 value_offset 0 - PCSCitationGeoKey: NAD_1983_HARN_StatePlane_Illinois_West_FIPS_1202_Feet|projection: Transverse
 Mercator
      key 3075 tiff_tag_location 0 count 1 value_offset 1 - ProjCoordTransGeoKey: CT_TransverseMercator
      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 5 - ProjLinearUnitSizeGeoKey: 0.3048006096
      key 3081 tiff_tag_location 34736 count 1 value_offset 0 - ProjNatOriginLatGeoKey: 36.66666667
      key 3082 tiff_tag_location 34736 count 1 value_offset 1 - ProjFalseEastingGeoKey: 2296583.333
      key 3083 tiff_tag_location 34736 count 1 value_offset 2 - ProjFalseNorthingGeoKey: 0
      key 3088 tiff_tag_location 34736 count 1 value_offset 3 - ProjCenterLongGeoKey: -90.16666667
      key 3092 tiff_tag_location 34736 count 1 value_offset 4 - ProjScaleAtNatOriginGeoKey: 0.9999411765
      key 4097 tiff_tag_location 34737 count 25 value_offset 151 - VerticalCitationGeoKey: NAVD88 - Geoid12B (Feet)
      key 4099 tiff_tag_location 0 count 1 value_offset 9003 - VerticalUnitsGeoKey: Linear_Foot_US_Survey
variable length header record 2 of 3:
  reserved             0
  user ID              'LASF_Projection'
  record ID            34736
  length after header  80
  description          'GeoTiff double parameters'
    GeoDoubleParamsTag (number of doubles 10)
      36.6667 2.29658e+006 0 -90.1667 0.999941 0.304801 6.37814e+006 298.257 0 0.0174533
variable length header record 3 of 3:
  reserved             0
  user ID              'LASF_Projection'
  record ID            34737
  length after header  290
  description          'GeoTiff ASCII parameters'
    GeoAsciiParamsTag (number of characters 290)
      NAD_1983_HARN_StatePlane_Illinois_West_FIPS_1202_Feet|projection: Transverse Mercator|PCS Name = NAD_1983_HARN_StatePlane_Illinois_West_FIPS_1202_Feet|NAV
D88 - Geoid12B (Feet)|GCS Name = GCS_North_American_1983_HARN|Datum = North_American_1983_HARN|Ellipsoid = GRS_1980|Primem = Greenwich|
the header is followed by 4 user-defined bytes
reporting minimum and maximum for all LAS point record entries ...
  X           568284000  573564000
  Y            73818000   79098000
  Z              596760     719240
  intensity           0      65535
  return_number       1          5
  number_of_returns   1          5
  edge_of_flight_line 0          1
  scan_direction_flag 0          0
  classification      1         15
  scan_angle_rank   -22         20
  user_data           0          1
  point_source_ID   140        142
  gps_time 101688614.763893 101690493.156573
WARNING: there is coordinate resolution fluff (x10) in XYZ
number of first returns:        12831281
number of intermediate returns: 152794
number of last returns:         12830952
number of single returns:       12214053
overview over number of returns of given pulse: 12214053 949868 395006 40987 1060 0 0
histogram of classification of points:
            1925  unclassified (1)
         7608322  ground (2)
         3089250  low vegetation (3)
          220295  medium vegetation (4)
          538780  high vegetation (5)
          164208  building (6)
             352  noise (7)
         1353327  wire conductor (14)
          624515  tower (15)
illinois_nad83_harn_sp_il_east.jpg
illinois_nad83_harn_sp_il_west.jpg

Kirk Waters - NOAA Federal

unread,
Feb 19, 2019, 7:23:46 AM2/19/19
to LAStools - efficient command line tools for LIDAR processing
Martin,
That's a good question. Unfortunately, you may actually need to further define the question. NAD83(HARN) is one of the NAD83 realizations. There are programs to transform between it and other NAD83 realizations on the NOAA NGS site. That part is reasonably well defined. The less defined part is which realization of WGS84 are we talking about? Different WGS84 realizations tend to match specific ITRF realizations. There are Helmert transformations published to go between NAD83 realizations and ITRF realizations. Below is a chunk of code I've got in my helmert class. The links in the comments may be what you're after.

/* 
* Initialize based on parameters we already know for certain transforms
*/
int helmert::init(char *from, char *to){
if(from == NULL || to == NULL) return 1;

/*
* The NAD83_86 transform values are from VDatum entries. Units are meters for
* tranlations and radians for rotations. No time dependencies.
* The ITRF2008/IGS08 to NAD83(2011) values are from https://www.ngs.noaa.gov/CORS/coords.shtml#Col2Exp on May 4, 2017.
* Units are meters for translations and milli arc seconds from rotations.
*
* NB: Comparison with HTDP indicates that HTDP does a transform relative to ITRF94. The ITRF94 to ITRF2008 values don't agree with the values from ITRF.
*/

if (
nameInit(from, to, "NAD83_86","WGS84_G1150", 0.9956, -1.9013, -0.52145, 1.256374e-7, 4.5701e-8, 5.6235e-8, 0.61504e-9,0) ||
nameInit(from, to, "NAD83_86","ITRF2000",0.9956, -1.9013, -0.52145, 1.256374e-7, 4.5701e-8, 5.6235e-8, 0.61504e-9,0) ||
nameInit(from, to, "NAD83_86","WGS84_G873",0.9910, -1.9072, -0.5129, 1.25033e-7, 4.6785e-8, 5.6529e-8, 0.0,0) ||
nameInit(from, to, "NAD83_86","ITRF96",0.9910, -1.9072, -0.5129, 1.25033e-7, 4.6785e-8, 5.6529e-8, 0.0,0) ||
nameInit(from, to, "NAD83_86","WGS84_G730",0.9830, -1.9092, -0.5049, 1.25033e-7, 4.6785e-8, 5.6529e-8, 8.0e-10,0) ||
nameInit(from, to, "NAD83_86","ITRF91",0.9830, -1.9092, -0.5049, 1.25033e-7, 4.6785e-8, 5.6529e-8, 8.0e-10,0) ||
nameInit(from, to, "NAD83_86","ITRF97",0.9956, -1.9013, -0.52145, 1.256374e-7, 4.5701e-8, 5.6235e-8, 0.61504e-9,0) ||
nameInit(from, to, "ITRF2008", "NAD83_MA11", 0.9080, -2.0161, -0.5653, 28.971, 10.420,  8.928, 1.10e-9, 0.0001, -0.0001, -0.0018, -0.020,  0.105, -0.347, 0.08e-9, 1997.0, 1) || // units in milli arc seconds for rotations
nameInit(from, to, "ITRF2008", "NAD83_2011", 0.99343, -1.90331, -0.52655, 25.91467,  9.42645, 11.59935, 1.71504e-9, 0.00079, -0.00060, -0.00134, 0.06667, -0.75744, -0.05133, -0.10201e-9, 1997.0, 1) ||
nameInit(from, to, "ITRF2008", "NAD83_PA11", 0.9080, -2.0161, -0.5653, 27.741, 13.469, 2.712, 1.10e-9, 0.0001, 0.0001, -0.0018, -0.384, 1.007, -2.186, 0.08e-9, 1997.0, 1) ||
nameInit(from, to, "ITRF2008", "ITRF2005", -2.0e-3, -0.9e-3, -4.7e-3, 0.0, 0.0, 0.0, 0.94e-9,  0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2000.0, 1) ||
nameInit(from, to, "ITRF2008", "ITRF2000", -1.9e-3, -1.7e-3, -10.5e-3, 0.00, 0.00, 0.00, 1.34e-9, 0.1e-3, 0.1e-3, -1.8e-3, 0.00, 0.00, 0.00, 0.08e-9, 2000.0,1) ||
nameInit(from, to, "ITRF2008", "ITRF97", 4.8e-3, 2.6e-3, -33.2e-3, 0.00, 0.00, 0.06, 2.92e-9, 0.1e-3, -0.5e-3, -3.2e-3, 0.00, 0.00, 0.02, 0.09e-9, 2000.0, 1) ||
nameInit(from, to, "ITRF2008", "ITRF96", 4.8e-3, 2.6e-3, -33.2e-3, 0.00, 0.00, 0.06, 2.92e-9, 0.1e-3, -0.5e-3, -3.2e-3, 0.00, 0.00, 0.02, 0.09e-9, 2000.0, 1) ||
nameInit(from, to, "ITRF2008", "ITRF94", 4.8e-3, 2.6e-3, -33.2e-3, 0.00, 0.00, 0.06, 2.92e-9, 0.1e-3, -0.5e-3, -3.2e-3, 0.00, 0.00, 0.02, 0.09e-9, 2000.0, 1) ||
nameInit(from, to, "ITRF2008", "NAD83_86", 0.99343 , -1.90331 , -0.52655 , -25.91467 , -9.42645 , -11.59935 , 1.71504e-9, 0.00079 , -0.00060 , -0.00134 , -0.06667 , 0.75744 , 0.05133 , -0.10201e-9, 1997.0, 1)
){

return 0;
}
//fprintf(stderr,"Didn't find a match for datums %s and %s\n",from,to);
return 1;
}

Kirk Waters, PhD, BJCP          | NOAA Office for Coastal Management
Applied Sciences Program      | 2234 South Hobson Ave
843-740-1227                          | Charleston, SC 29405    



Evon Silvia

unread,
Feb 19, 2019, 6:13:59 PM2/19/19
to last...@googlegroups.com
Officially you're supposed to use NGS's GEOCON tool for the HARN transformations, which is now part of the beta NGS NCAT tool (https://www.ngs.noaa.gov/NCAT/).

Unfortunately, it's not usually that simple because as Kirk points out there are multiple WGS84 realizations and often the HARN data isn't very well defined. Usually for web visualization it's "close enough" (+/- a couple meters) to just de-project the data and hold NAD83 and WGS84 equivalent.

If you want to go more precise than that it's a can of worms and depends on how the project was controlled. I cringe every time a client requests NAD83(HARN) for this reason. At least NAD83(2011) was well defined and a national system, but HARN was defined independently for every state, sometimes twice as in in the Pacific Northwest where I live...

Evon
Reply all
Reply to author
Forward
0 new messages