Multispectral LiDAR (light detection and ranging) data have been initially used for land cover classification. However, there are still high classification uncertainties, especially in urban areas, where objects are often mixed and confounded. This study investigated the efficiency of combining advanced statistical methods and LiDAR metrics derived from multispectral LiDAR data for improving land cover classification accuracy in urban areas. The study area is located in Oshawa, Ontario, Canada, on the Lake Ontario shoreline. Multispectral Optech Titan LiDAR data over the study area were acquired on 3 September 2014 in a single strip of 3 km2. Using the channels at 1,550 nm (C1), 1,064 nm (C2) and 532 nm (C3), LiDAR intensity data, normalized digital surface model (nDSM), pseudo normalized difference vegetation index (PseudoNDVI), morphological profiles (MP), and a novel hierarchical morphological profiles (HMP) were derived and used as features for the classification. A support vector machine classifier with a radial basis function (RBF) kernel was applied in the classification stage, where the optimal parameters for the classifier were selected by a grid search procedure. The combination of intensity, pseudoNDVI, nDSM and HMP resulted in the best land cover classification, with an overall accuracy of 93.28%.