[GRASS-user] Need help importong LiDAR file with r.in.gdal

rshepard · May 1, 2023, 9:01pm

I'm trying to import a LiDAR data set. I've used EPSG:6559 while the data
source manager tells me it imports into ARC/gis using EPSG:6557. I've tried
both and grass doesn't work for me regardless of which code I use.

I create a location:
grass -c EPSG:6557 /data1/grassdata/Metro_LiDAR

Then import the data set (N.B. a copy of the data source tarball is on the
fileconvoy.com web site; see below):
r.in.gdal in=Bare_Earth/be45122f8/hdr.adf out=topo_2010 mem=2000 --o

GRASS returns:
ERROR: Projection of dataset does not appear to match current location.

        Location PROJ_INFO is:
        name: NAD83(2011) / Oregon GIC Lambert (ft)
        datum: nad83_2011
        ellps: grs80
        proj: lcc
        lat_0: 41.75
        lon_0: -120.5
        lat_1: 43
        lat_2: 45.5
        x_0: 399999.9999984
        y_0: 0
        no_defs: defined

        Dataset PROJ_INFO is:
        name: unnamed
        ellps: grs80
        proj: lcc
        lat_0: 41.75
        lon_0: -120.5
        lat_1: 43
        lat_2: 45.5
        x_0: 400000
        y_0: 0
        no_defs: defined

Difference in: datum

        In case of no significant differences in the projection definitions,
        use the -o flag to ignore them and use current location definition.
        Consider generating a new location from the input dataset using the
        'location' parameter.

I added the -o option to r.in.gdal and this is the result:
r.in.gdal -o in=$HOME/projects/oregon/northside-rock/data/LDQ-45122F8/2014-OLC-Metro/Bare_Earth/be45122f8/hdr.adf out=topo_2010 mem=2000 --o
Over-riding projection check
Importing raster map <topo_2010>...
ERROR 1: Block is corrupt, block size was 21532, but expected to be 16384.
ERROR 1: /home/rshepard/projects/oregon/northside-rock/data/LDQ-45122F8/2014-OLC-Metro/Bare_Earth/be45122f8, band 1: IReadBlock failed at X offset 15, Y offset 308: Block is corrupt, block size was 21532, but expected to be 16384.
ERROR: Unable to read row 4928

If the file loads into ARC/gis why not into grass?

I do not see why this isn't working here. I've uploaded the Bare_Earth
directory to
http://www.fileconvoy.com/dfl.php?id=g4e3d18b830b60bef10004916731ceed9a68e0d28c3.
It's a 474M file.

I need to learn why I have problems importing the data using either EPSG
code, ignoring the data project by override or not.

Regards,

Rich

neteler · May 1, 2023, 9:21pm

Rich,

On Mon, May 1, 2023 at 11:01 PM Rich Shepard <rshepard@appl-ecosys.com> wrote:

I'm trying to import a LiDAR data set. I've used EPSG:6559 while the data
source manager tells me it imports into ARC/gis using EPSG:6557. I've tried
both and grass doesn't work for me regardless of which code I use.

[...]

r.in.gdal -o in=$HOME/projects/oregon/northside-rock/data/LDQ-45122F8/2014-OLC-Metro/Bare_Earth/be45122f8/hdr.adf out=topo_2010 mem=2000 --o
Over-riding projection check
Importing raster map <topo_2010>...
ERROR 1: Block is corrupt, block size was 21532, but expected to be 16384.
ERROR 1: /home/rshepard/projects/oregon/northside-rock/data/LDQ-45122F8/2014-OLC-Metro/Bare_Earth/be45122f8, band 1: IReadBlock failed at X offset 15, Y offset 308: Block is corrupt, block size was 21532, but expected to be 16384.
ERROR: Unable to read row 4928

What does it show in GDAL?

gdalinfo $HOME/projects/oregon/northside-rock/data/LDQ-45122F8/2014-OLC-Metro/Bare_Earth/be45122f8/hdr.adf

The point is that GDAL needs to be able to read it (since r.in.gdal
depends on it).

Markus

rshepard · May 1, 2023, 9:33pm

On Mon, 1 May 2023, Markus Neteler wrote:

What does it show in GDAL?

Markus,

I forgot about gdalinfo and the file to which it's applied.

For this data set:
$ gdalinfo hdr.adf
Driver: AIG/Arc/Info Binary Grid
Files: .
        ./prj.adf
        ./w001001.adf
        ./w001001x.adf
        ./sta.adf
        ./metadata.xml
        ./dblbnd.adf
        ./hdr.adf
Size is 11097, 15505
Coordinate System is:
PROJCRS["unnamed",
     BASEGEOGCRS["Unknown datum based upon the GRS 1980 ellipsoid",
         DATUM["Not specified (based on GRS 1980 ellipsoid)",
             ELLIPSOID["GRS 1980",6378137,298.257222101,
                 LENGTHUNIT["metre",1]]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433]],
         ID["EPSG",4019]],
     CONVERSION["unknown",
         METHOD["Lambert Conic Conformal (2SP)",
             ID["EPSG",9802]],
         PARAMETER["Latitude of false origin",41.75,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8821]],
         PARAMETER["Longitude of false origin",-120.5,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8822]],
         PARAMETER["Latitude of 1st standard parallel",43,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8823]],
         PARAMETER["Latitude of 2nd standard parallel",45.5,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8824]],
         PARAMETER["Easting at false origin",1312335.95800525,
             LENGTHUNIT["user-defined",0.3048],
             ID["EPSG",8826]],
         PARAMETER["Northing at false origin",0,
             LENGTHUNIT["user-defined",0.3048],
             ID["EPSG",8827]]],
     CS[Cartesian,2],
         AXIS["(E)",east,
             ORDER[1],
             LENGTHUNIT["user-defined",0.3048]],
         AXIS["(N)",north,
             ORDER[2],
             LENGTHUNIT["user-defined",0.3048]]]
Data axis to CRS axis mapping: 1,2
Origin = (672763.500000000000000,1467820.500000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Corner Coordinates:
Upper Left ( 672763.500, 1467820.500) (123d 0'19.54"W, 45d44'59.57"N)
Lower Left ( 672763.500, 1421305.500) (122d59'59.58"W, 45d37'20.60"N)
Upper Right ( 706054.500, 1467820.500) (122d52'30.33"W, 45d45' 9.33"N)
Lower Right ( 706054.500, 1421305.500) (122d52'11.40"W, 45d37'30.34"N)
Center ( 689409.000, 1444563.000) (122d56'15.21"W, 45d41'15.03"N)
Band 1 Block=256x16 Type=Float32, ColorInterp=Undefined
   Min=11.370 Max=1163.480
   NoData Value=-3.4028235e+38

Yet ARC/gis apparently has no issues with it.

Regards,

Rich

rshepard · May 2, 2023, 4:39pm

On Mon, 1 May 2023, Rich Shepard wrote:

Coordinate System is:
PROJCRS["unnamed",
   BASEGEOGCRS["Unknown datum based upon the GRS 1980 ellipsoid",
       DATUM["Not specified (based on GRS 1980 ellipsoid)",
           ELLIPSOID["GRS 1980",6378137,298.257222101,
               LENGTHUNIT["metre",1]]],
       PRIMEM["Greenwich",0,
           ANGLEUNIT["degree",0.0174532925199433]],
       ID["EPSG",4019]],
   CONVERSION["unknown",
       METHOD["Lambert Conic Conformal (2SP)",
           ID["EPSG",9802]],
       PARAMETER["Latitude of false origin",41.75,
           ANGLEUNIT["degree",0.0174532925199433],
           ID["EPSG",8821]],
       PARAMETER["Longitude of false origin",-120.5,
           ANGLEUNIT["degree",0.0174532925199433],
           ID["EPSG",8822]],
       PARAMETER["Latitude of 1st standard parallel",43,
           ANGLEUNIT["degree",0.0174532925199433],
           ID["EPSG",8823]],
       PARAMETER["Latitude of 2nd standard parallel",45.5,
           ANGLEUNIT["degree",0.0174532925199433],
           ID["EPSG",8824]],
       PARAMETER["Easting at false origin",1312335.95800525,
           LENGTHUNIT["user-defined",0.3048],
           ID["EPSG",8826]],
       PARAMETER["Northing at false origin",0,
           LENGTHUNIT["user-defined",0.3048],
           ID["EPSG",8827]]],
   CS[Cartesian,2],
       AXIS["(E)",east,
           ORDER[1],
           LENGTHUNIT["user-defined",0.3048]],
       AXIS["(N)",north,
           ORDER[2],
           LENGTHUNIT["user-defined",0.3048]]]
Data axis to CRS axis mapping: 1,2
Origin = (672763.500000000000000,1467820.500000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Corner Coordinates:
Upper Left ( 672763.500, 1467820.500) (123d 0'19.54"W, 45d44'59.57"N)
Lower Left ( 672763.500, 1421305.500) (122d59'59.58"W, 45d37'20.60"N)
Upper Right ( 706054.500, 1467820.500) (122d52'30.33"W, 45d45' 9.33"N)
Lower Right ( 706054.500, 1421305.500) (122d52'11.40"W, 45d37'30.34"N)
Center ( 689409.000, 1444563.000) (122d56'15.21"W, 45d41'15.03"N)

If these are valid values why might grass have problems importing the map?

Regards,

Rich