This work introduces range-based GP maps, which directly represent terrain by modeling the range from a LiDAR sensor as a Gaussian process (GP) in spherical space. Such a model aligns the predicted uncertainty from the GP regression with the uncertainty in the underlying sensor observations. Experimental evaluation on simulated natural terrain indicates that local range-based GP maps perform comparably to elevation-based methods when predicting terrain height, with the former producing more stable parameters and providing a better uncertainty representation. An aggregation method is proposed using the pose as an additional input to the GP. Unlike their elevation-based counterparts, range-based GP maps are capable of modeling overhangs and vertical obstacles with ease, demonstrated with examples of maps built on real-world data from a fully 3D subterranean environment.