EPSG code 1256 not found in 'vertcs.csv' file

[Spencer Floyd]

Hello,

 
When I run the following command:

blast2dem -i *.laz -cores 7 -keep_class 2 -step 1 -use_tile_bb - obil

I'm getting the following error:

Below is the WKT OGC CRS taken from lasinfo when run on one of the .laz files. I have bolded the epsg in question:

WKT OGC COORDINATE SYSTEM:
    COMPD_CS["NAD83(CSRS)v7 / UTM zone 10N + CGVD2013a(2010) height",PROJCS["NAD83(CSRS)v7 / UTM zone 10N",GEOGCS["NAD83(CSRS)v7",DATUM["North_American_Datum_of_1983_CSRS_version_7",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1198"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","8255"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","22710"]],VERT_CS["CGVD2013a(2010) height",VERT_DATUM["Canadian Geodetic Vertical Datum of 2013 (CGG2013a) epoch 2010",2005,AUTHORITY["EPSG","1256"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","9245"]]

This seems to be the DATUM_CODE when I look for it in the vertcs.csv

The real EPSG is 9245.

It seems to exist in the vertcs.csv, so why is this error happening?

I use -load_ogc_wkt to bring back the compound WTK projection using a .txt file which contains the correct WKT CRS string. 

Thanks,

Spencer

[Jochen Rapidlasso]

Hi Spencer,

after doing some tests with this wkt:

It seems the wkt parser just struggles with the total string length of almost 1000 characters. If you separate the wkt into multiple lines and add the last missing closing bracket,
blast2dem will take it without any errors.

This one worked:

COMPD_CS["NAD83(CSRS)v7 / UTM zone 10N + CGVD2013a(2010) height",PROJCS["NAD83(CSRS)v7 / UTM zone 10N",GEOGCS["NAD83(CSRS)v7",DATUM["North_American_Datum_of_1983_CSRS_version_7",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1198"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","8255"]],PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","22710"]],

VERT_CS["CGVD2013a(2010) height",VERT_DATUM["Canadian Geodetic Vertical Datum of 2013 (CGG2013a) epoch 2010",2005,AUTHORITY["EPSG","1256"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","9245"]]]

Best,

Jochen @rapidlasso

[Spencer Floyd]

Hey Jochen,

Thanks for the response. I assume you're referring to separating the wkt into multiple lines within the text file I'm using with -load_ogc_wkt in las2las.

I tried this and when I run lasinfo, it seems that only the first line gets used to create the new WKT CRS:


You're right that it can be run through Blast2DEM now, but when you run gdalinfo on the resulting raster, there is no CRS information:

[Jochen Rapidlasso]

Hi Spencer,

you are right, that was a quick shot, and a wrong one.

After debugging into it:

We look for the first column in vertcs.csv for the crs epsg (1st column), out of your WKT we fetch the datum.

To solve the problem in your situation it is probably the best just to give this values straight as argument:

    blast2dem -i *.laz -cores 7 -keep_class 2 -step 1 -use_tile_bb - obil -epsg 22710 -vertical_epsg 9245

should do the job.

Cheers,

Jochen @rapidlasso

[Spencer Floyd]

Hey Jochen,

I tried using the -epsg and -vertical_epsg parameters.

This is the output I get in laslook:

When I look at the the gdalinfo output run on the resulting .tif file we have the horizontal datum but no vertical datum:

I think it may due to the fact that these two particular datums don't have a compound projection epsg code (that I've been able to find). 

I'm certainly no expert in datums and projections, I'm just trying to make both these horizontal and vertical datums work.

If anyone has any advice, I'd appreciate it.

Thanks,

Spencer

[Michael Perdue]

Hi Spencer,

I don’t think any lastool program embes the vertical datum in outputs unless (as you guessed) the EPSG code defines a compound datum.

This is probably the better default behavior as most of the tiff outputs are not elevation data. lasoverlap, lascanopy and lasgrid all output files that are not elevation data. Lasgrid might be lowest point, highest point or average height, but it might also be any one of the other fields or metrics such as intensity or density. These are data that have no vertical component and hence, no vertical datum (IMHO). Even las2dem can produce outputs that should not have a vertical datum applied.

To implement it means that LAStools needs to be smart and know when the output should and should not include the vertical datum. Maybe Rapidlasso will add that in the future, but in the meantime you can apply the proper SRS information to your output as necessary through gdal_edit.py. The edit updates the header inline so the task is fast. Even tens of thousands of tiles should finish in a few minutes.

michaelperdue@Hactar /mnt/e/tmp/density_geotiffs

>gdal_edit.py -a_srs "EPSG:22710+9245" 7005_10100.tif

michaelperdue@Hactar /mnt/e/tmp/density_geotiffs

>gdalinfo 7005_10100.tif

Driver: GTiff/GeoTIFF

Files: 7005_10100.tif

Size is 250, 232

Coordinate System is:

COMPOUNDCRS["NAD83(CSRS)v7 / UTM zone 10N + CGVD2013a(2010) height",

PROJCRS["NAD83(CSRS)v7 / UTM zone 10N",

BASEGEOGCRS["NAD83(CSRS)v7",

DATUM["North American Datum of 1983 (CSRS) version 7",

ELLIPSOID["GRS 1980",6378137,298.257222101,

LENGTHUNIT["metre",1]]],

PRIMEM["Greenwich",0,

ANGLEUNIT["degree",0.0174532925199433]],

ID["EPSG",8255]],

CONVERSION["UTM zone 10N",

METHOD["Transverse Mercator",

ID["EPSG",9807]],

PARAMETER["Latitude of natural origin",0,

ANGLEUNIT["degree",0.0174532925199433],

ID["EPSG",8801]],

PARAMETER["Longitude of natural origin",-123,

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["Canada between 126°W and 120°W, onshore and offshore south of 84°N - British Columbia, Northwest Territories, Yukon."],

BBOX[48.13,-126,81.8,-120]],

ID["EPSG",22710]],

VERTCRS["CGVD2013a(2010) height",

VDATUM["Canadian Geodetic Vertical Datum of 2013 (CGG2013a) epoch 2010",

ANCHOREPOCH[2010]],

CS[vertical,1],

AXIS["gravity-related height (H)",up,

LENGTHUNIT["metre",1]],

USAGE[

SCOPE["Geodesy, topographic mapping."],

AREA["Canada - onshore and offshore - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon."],

BBOX[38.21,-141.01,86.46,-40.73]],

ID["EPSG",9245]]]

Data axis to CRS axis mapping: 1,2,3

Origin = (700500.000000000000000,-1009536.000000000000000)

Pixel Size = (2.000000000000000,-2.000000000000000)

Metadata:

AREA_OR_POINT=Area

TIFFTAG_ARTIST=created by LAStools (c) in...@rapidlasso.de

TIFFTAG_SOFTWARE=LAStools is a product of rapidlasso GmbH, Germany

Image Structure Metadata:

COMPRESSION=LZW

INTERLEAVE=BAND

Corner Coordinates:

Upper Left ( 700500.000,-1009536.000) (121d10'31.78"W, 9d 7'41.64"S)

Lower Left ( 700500.000,-1010000.000) (121d10'31.70"W, 9d 7'56.74"S)

Upper Right ( 701000.000,-1009536.000) (121d10'15.40"W, 9d 7'41.56"S)

Lower Right ( 701000.000,-1010000.000) (121d10'15.33"W, 9d 7'56.66"S)

Center ( 700750.000,-1009768.000) (121d10'23.55"W, 9d 7'49.15"S)

Band 1 Block=250x20 Type=Float32, ColorInterp=Gray

NoData Value=0

Unit Type: metre

Note: This is an example where the file should not have a vertical datum. The values in these files are density, but were used for illustration.

--
Download LAStools at
https://rapidlasso.de
Manage your settings at
https://groups.google.com/g/lastools/membership
---
You received this message because you are subscribed to the Google Groups "LAStools - efficient tools for LiDAR processing" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lastools+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lastools/2f714313-5a11-4f07-8f2d-498a05456153n%40googlegroups.com.

[Spencer Floyd]

Michael,

Thanks so much for posting this. That seems like the way to go at this point.

Spencer