A robust and efficient method of calculating a cartesian grid of heights or roughnesses from contour line maps is developed. The purpose of the grids is to serve as input for atmospheric flow solvers such as WAsP Engineering or EllipSys3D. The method builds on Delaunay triangulation constrained to include all contour segments in the triangulation. It is furthermore refined to avoid spurious flat areas produced by the Delaunay triangulation. Robust ways to extrapolate beyond the convex hull of the map points are provided.