processing lidar with liblas and gdal

kivalina-lidar-process.sh

# depends on liblas built with laszip and ogr support

mkdir -p laz
cd laz
# download from noaa
if [ ! -f 20040805A_ak302_2_000037_p13.laz ]; then
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000037_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000038_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000039_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000040_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000041_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000042_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000043_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000044_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000045_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000046_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000047_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000049_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000050_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000051_p13.laz
  wget http://www.csc.noaa.gov/htdata/lidar1_z/geoid12a/data/23/20040805A_ak302_2_000052_p13.laz
fi
cd ../

mkdir -p las
# uncompress laz to las
for i in $(ls laz/*laz); do
  las2las -i $i -o $i.las
  mv $i.las las/
done

mkdir -p shp 
# convert las to shp
for i in $(ls las/*las); do
  las2ogr -i $i -o $i.geo.shp
  mv $i.geo.* shp/
done
 
# set up projections as variables
export SOURCE_EPSG="EPSG:4269"
export UTM_EPSG="EPSG:3470" # NAD83 Alaska zone 3
export MERC_EPSG="EPSG:3857"
export ALBERS_EPSG="EPSG:3467"
 
# merge all shapefiles into one
rm merged.*
ogr2ogr merged.shp shp/20040805A_ak302_2_000045_p13.laz.las.geo.shp
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000037_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000038_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000039_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000040_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000041_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000042_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000043_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000044_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000045_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000046_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000047_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000049_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000050_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000051_p13.laz.las.geo.shp -update -append
ogr2ogr merged.shp -nln merged shp/20040805A_ak302_2_000052_p13.laz.las.geo.shp -update -append
 
# import into postgis for an easy way to throw away data
createdb -T template_postgis points -h localhost
time shp2pgsql -s 4269 merged.shp merged | psql points -h localhost -qq
# real  6m19.752s

# find anomalies
psql -h localhost points -c "select ST_Z(geom) from merged  order by ST_Z(geom) desc limit 15;"
: '
   st_z   
----------
 135.1136
 132.0859
  82.5735
  77.4617
  65.2517
  57.7734
  53.5563
  39.4419
  36.0878
  36.0878
  23.9503
  23.9503
  23.7603
  23.7603
  23.7304
(15 rows)
'

# then dump back out while clipping values to what look like ground and throwing away high values that are clearly anomalies
pgsql2shp -f merged_over_2_under_22.shp points 'SELECT gid, geom, ST_Z(geom) as elev, intensity from merged where ST_Z(geom) >= 2.1 and ST_Z(geom) <= 21'
 
# build spatial index
ogrinfo merged_over_2_under_22.shp -sql "CREATE SPATIAL INDEX ON merged_over_2_under_22"
 
# make overall bounds shapefile
ogrtindex merged_index.shp merged_over_2_under_22.shp
 
# clip to kivalina proper
time ogr2ogr -clipsrc -164.56471 67.7185 -164.51915 67.7348 merged-points-wgs84.shp merged_over_2_under_22.shp
 
# NOW manually edit in QGIS to remove a few odd located points
 
# create a dem using gdal_grid
# todo digest http://trac.osgeo.org/gdal/ticket/2411
SHP="merged-points-wgs84.shp"
LAYER="merged-points-wgs84"
SIZE="1000"
FIELD="ELEV"
TYPE="dem"
OUT="kiv-${TYPE}-wgs84-tmp-${SIZE}.tif"
OUT2="kiv-${TYPE}-wgs84.tif"
time gdal_grid -a nearest -zfield ${FIELD} -a_srs $SOURCE_EPSG -outsize ${SIZE} ${SIZE} ${SHP} -l ${LAYER} ${OUT} --config GDAL_NUM_THREADS ALL_CPUS
 
# fix broken geotransform 
gdalwarp ${OUT} ${OUT2}
 
# generate contours
time gdal_contour -a z $OUT2 kiv-contour-.5-wgs84.shp -i .5