#!/bin/bash
# Konsta Happonen, 2015

# Run this script from within GRASS 7 compiled with libLAS.
# It creates variables of forest structure from its input files
# The input files are supposed to be georeferenced .las files.
#
# You need to have imported the utm5 mapsheet division into GRASS's PERMANENT map with the name utm5. 
#
# the directory ASCIINAME must exist
#
# You will also need the DEM that corresponds to the utm5-sheet of your .las file.
#
# This script will not work out of the box. You need to know what you're doing to modify it
# and make it work.

LASFILE=$1
BASEFILE=${LASFILE##*/}
BASEFILE=${BASEFILE%.las}
DEMPATH="/mnt/grassdata/EPSG:3067/PERMANENT/laserkeilaus/demit/"
DEMBASEFILE=${BASEFILE%?}
DEM=${DEMPATH}${DEMBASEFILE}".asc"
ASCIINAME="/mnt/grassdata/EPSG:3067/PERMANENT/laserkeilaus/ascii/${BASEFILE}.txt" 
POINT_MAP="${BASEFILE}_points"
INT_RASTER="${BASEFILE}_intensity_2m"
MAXH_RASTER="${BASEFILE}_max_height_6m"
DEM_RASTER="${BASEFILE}_DEM_2m"
NONGROUND_RASTER="${BASEFILE}_nonground_n_6m"
N_RASTER="${BASEFILE}_n_6m"
COVER_RASTER="${BASEFILE}_canopycover_6m"

echo "LASFILE: $LASFILE"
echo "BASEFILE: $BASEFILE"
echo "DEMPATH: $DEMPATH"
echo "DEMBASEFILE: $DEMBASEFILE"
echo "DEM: $DEM"
echo "ASCIINAME: $ASCIINAME"
echo "POINT_MAP: $POINT_MAP"
echo "INT_RASTER: $INT_RASTER"
echo "MAXH_RASTER: $MAXH_RASTER"
echo "DEM_RASTER: $DEM_RASTER"
echo "NONGROUND_RASTER: $NONGROUND_RASTER"
echo "N_RASTER: $N_RASTER"
echo "COVER_RASTER: $COVER_RASTER"

# Check that the .las file and DEM exist.
if [ -f $LASFILE ];
then
    echo "$LASFILE exists"
else
    echo "$LASFILE doesn't exist."
    echo "Exiting..."
    exit 1
fi

if [ -f $DEM ];
then
    echo "$DEM exists"
else
    echo "$DEM doesn't exist."
    echo "Exiting..."
    exit 1
fi

# Transform the .las file into a point vector feature, because it's easier to handle in GRASS.
echo "Transforming .las-file into point vector form..."
v.in.lidar -o --overwrite input=$LASFILE output=$POINT_MAP class_filter=1,2,3
echo "Done."

# Create new columns into the point vector's datasheet, DEM (ground elevation) and normalized DSM ("canopy height")
echo "Adding columns into the point feature's datasheet for DEM and NDSM..."
v.db.addcolumn map=$POINT_MAP columns="dem_value double precision,ndsm_value double precision"
echo "Done."

# Set GRASS region based on the .las file.
echo "Setting region..."
v.extract --overwrite input=utm5@keskisuomi output=$BASEFILE where="LEHTITUNNU = \"$BASEFILE\""
g.region vector=$BASEFILE res=2 -ap
echo "Done."

# Interpolating intensity into a 2 m resolution raster with a smoothing function based on inverse distance
echo "Creating intensity raster..."
v.surf.idw --overwrite input=$POINT_MAP column=intensity output=$INT_RASTER npoints=12 power=12
echo "Done."

# Import DEM into GRASS. CRS should be ETRS-TM35FIN.
echo "Importing DEM..."
r.in.gdal --overwrite -o input=$DEM output=$DEM_RASTER
echo "Done"

# Fetch elevations from DEM to the point feature's datasheet
echo "Sampling elevations to the point feature's datasheet..."
v.what.rast map=$POINT_MAP column=dem_value raster=$DEM_RASTER
echo "Done."

# Calculate NDSM-values for the points
echo "Calculating NDSM..."
v.db.update map=$POINT_MAP column=ndsm_value query_column="z_coord-dem_value"
echo "Done."

# This is what it looks like right now:
# cat|x_coord|y_coord|z_coord|intensity|return|n_returns|scan_dir|edge|cl_type|class|gps_time|angle|src_id|usr_data|red|green|blue|GRASSRGB|dem_value|ndsm_value
# 1|377499.94|6897012.98|106.37|11|1|1|0|0|Synthetic|Water|51464675.258683|-9|7|1|0|0|0|000:000:000|106.394|-0.0240000000000009


# To use r.in.xyz, the feature must be exported as text :/
echo "Exporting point feature as text for rasterzing..."
v.out.ascii --overwrite input=$POINT_MAP type=point output=$ASCIINAME columns="ndsm_value"
echo "Done."

# Lasketaan 6 m resoluutiolla maksimikorkeusrasteri normalisoidusta pintamallista (ndsm_value)
# Monesko sarake ASCII-tiedostossa on NDSM?
# 377499.94|6897012.98|106.37|1|-0.02400000
# 5. sarake.
# Se laitetaan z-arvoksi.

echo "Setting resolution for maximum height raster"
g.region vector=$BASEFILE res=6 -ap

echo "Calculating maximum height raster..."
r.in.xyz --overwrite input=$ASCIINAME output=$MAXH_RASTER method=max type=FCELL z=5
echo "Done"

echo "Creating a raster from lidar returns classified as nonground..."
r.in.lidar --overwrite -o input=$LASFILE output=$NONGROUND_RASTER method=n resolution=6 type=CELL class_filter=1,3
echo "Done."

echo "Creating a raster from all returns..."
r.in.lidar --overwrite -o input=$LASFILE output=$N_RASTER method=n resolution=6 type=CELL class_filter=1,2,3
echo "Done."

echo "Calculating canopy cover index raster..."
r.mapcalc --overwrite expression="$COVER_RASTER = float(${NONGROUND_RASTER})/${N_RASTER}"
echo "Done."

echo "All done."

exit 0
