Skip to content

Instantly share code, notes, and snippets.

@springmeyer
Last active October 17, 2022 13:48
Show Gist options
  • Save springmeyer/6016995 to your computer and use it in GitHub Desktop.
Save springmeyer/6016995 to your computer and use it in GitHub Desktop.
processing lidar with liblas and gdal
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment