We consider the melting of the vortex solid in highly anisotropic layered superconductors with a small concentration of random columnar pinning centers. Using large-scale numerical minimization of a free-energy functional, we find that melting of the low-temperature, nearly crystalline vortex solid (Bragg glass) into a vortex liquid occurs in two steps as the temperature increases: the Bragg glass and liquid phases are separated by an intermediate Bose glass phase. A suitably defined local melting temperature exhibits spatial variation similar to that observed in experiments.