diff --git a/BGL/include/CGAL/boost/graph/IO/PLY.h b/BGL/include/CGAL/boost/graph/IO/PLY.h
index f14ad0173af..6cf64f5dd1f 100644
--- a/BGL/include/CGAL/boost/graph/IO/PLY.h
+++ b/BGL/include/CGAL/boost/graph/IO/PLY.h
@@ -335,7 +335,7 @@ bool write_PLY(std::ostream& os,
 
   typedef typename CGAL::GetInitializedVertexIndexMap<Graph, CGAL_BGL_NP_CLASS>::const_type VIMap;
   typedef typename GetVertexPointMap<Graph, CGAL_BGL_NP_CLASS>::const_type                  Vpm;
-  typedef typename boost::property_traits<Vpm>::reference                                   Point_3;
+  typedef typename boost::property_traits<Vpm>::value_type                                  Point_3;
   typedef CGAL::IO::Color                                                                   Color;
   typedef typename internal_np::Lookup_named_param_def<
                      internal_np::vertex_color_map_t,
@@ -407,7 +407,7 @@ bool write_PLY(std::ostream& os,
 
   for(vertex_descriptor vd : vertices(g))
   {
-    Point_3 p = get(vpm, vd);
+    const Point_3& p = get(vpm, vd);
     internal::output_properties(os, &p, make_ply_point_writer (CGAL::Identity_property_map<Point_3>()));
     if(has_vcolor)
     {
diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h
index ff690075698..0d8ca71e7d6 100644
--- a/BGL/include/CGAL/boost/graph/named_params_helper.h
+++ b/BGL/include/CGAL/boost/graph/named_params_helper.h
@@ -309,23 +309,26 @@ CGAL_DEF_GET_INITIALIZED_INDEX_MAP(face, typename boost::graph_traits<Graph>::fa
 
   template<typename PointRange,
            typename NamedParameters = Named_function_parameters<bool, internal_np::all_default_t>,
-           bool has_nested_iterator = internal::Has_nested_type_iterator<PointRange>::value>
+           bool has_nested_iterator = internal::Has_nested_type_iterator<PointRange>::value,
+           typename NP_TAG = internal_np::point_t
+  >
   class GetPointMap
   {
     typedef typename std::iterator_traits<typename PointRange::iterator>::value_type Point;
     typedef typename CGAL::Identity_property_map<Point> DefaultPMap;
+    typedef typename CGAL::Identity_property_map<const Point> DefaultConstPMap;
 
   public:
     typedef typename internal_np::Lookup_named_param_def<
-    internal_np::point_t,
+    NP_TAG,
     NamedParameters,
     DefaultPMap
     > ::type  type;
 
     typedef typename internal_np::Lookup_named_param_def<
-    internal_np::point_t,
+    NP_TAG,
     NamedParameters,
-    DefaultPMap
+    DefaultConstPMap
     > ::type  const_type;
   };
 
@@ -336,7 +339,7 @@ CGAL_DEF_GET_INITIALIZED_INDEX_MAP(face, typename boost::graph_traits<Graph>::fa
     struct Dummy_point{};
   public:
     typedef typename CGAL::Identity_property_map<Dummy_point> type;
-    typedef typename CGAL::Identity_property_map<Dummy_point> const_type;
+    typedef typename CGAL::Identity_property_map<const Dummy_point> const_type;
   };
 
   namespace Point_set_processing_3 {
@@ -375,26 +378,6 @@ CGAL_DEF_GET_INITIALIZED_INDEX_MAP(face, typename boost::graph_traits<Graph>::fa
         >::Kernel::FT type;
     };
 
-    template<typename PointRange, typename NamedParameters>
-    class GetQueryPointMap
-    {
-      typedef typename std::iterator_traits<typename PointRange::iterator>::value_type Point;
-      typedef typename CGAL::Identity_property_map<Point> DefaultPMap;
-
-    public:
-      typedef typename internal_np::Lookup_named_param_def<
-      internal_np::query_point_t,
-      NamedParameters,
-      DefaultPMap
-      > ::type  type;
-
-      typedef typename internal_np::Lookup_named_param_def<
-      internal_np::query_point_t,
-      NamedParameters,
-      DefaultPMap
-      > ::type  const_type;
-    };
-
     template<typename PointRange, typename NamedParameters>
     class GetK
     {
@@ -440,6 +423,7 @@ CGAL_DEF_GET_INITIALIZED_INDEX_MAP(face, typename boost::graph_traits<Graph>::fa
     {
       typedef typename PlaneRange::iterator::value_type Plane;
       typedef typename CGAL::Identity_property_map<Plane> DefaultPMap;
+      typedef typename CGAL::Identity_property_map<const Plane> DefaultConstPMap;
 
     public:
       typedef typename internal_np::Lookup_named_param_def<
@@ -451,7 +435,7 @@ CGAL_DEF_GET_INITIALIZED_INDEX_MAP(face, typename boost::graph_traits<Graph>::fa
       typedef typename internal_np::Lookup_named_param_def<
       internal_np::plane_t,
       NamedParameters,
-      DefaultPMap
+      DefaultConstPMap
       > ::type  const_type;
     };
 
diff --git a/Classification/include/CGAL/Classification/Point_set_neighborhood.h b/Classification/include/CGAL/Classification/Point_set_neighborhood.h
index bc1b2b64ef4..c204a337454 100644
--- a/Classification/include/CGAL/Classification/Point_set_neighborhood.h
+++ b/Classification/include/CGAL/Classification/Point_set_neighborhood.h
@@ -73,7 +73,9 @@ class Point_set_neighborhood
     My_point_property_map (const PointRange *input, PointMap point_map)
       : input (input), point_map (point_map) { }
 
-    friend reference get (const My_point_property_map& ppmap, key_type i)
+    // we did not put `reference` here on purpose as the recommanded default
+    // is `Identity_property_map<Point_3>` and not `Identity_property_map<const Point_3>`
+    friend decltype(auto) get (const My_point_property_map& ppmap, key_type i)
     { return get(ppmap.point_map, *(ppmap.input->begin()+std::size_t(i))); }
   };
 
diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h
index fc243904199..f6eb57c4217 100644
--- a/Point_set_processing_3/include/CGAL/estimate_scale.h
+++ b/Point_set_processing_3/include/CGAL/estimate_scale.h
@@ -250,7 +250,7 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_2
 
     friend inline value_type get (const Pmap_to_3d& pmap, const key_type& p)
     {
-      typename boost::property_traits<PointMap>::reference p2 = get(pmap.point_map, p);
+      const typename boost::property_traits<PointMap>::value_type& p2 = get(pmap.point_map, p);
       return value_type (p2.x(), p2.y(), 0.);
     }
   };
@@ -383,7 +383,7 @@ public:
     FT nb = 0.;
     std::size_t index = 0;
 
-    typename boost::property_traits<PointMap>::reference
+    const typename boost::property_traits<PointMap>::value_type&
       pquery = get(point_map, *query);
     for (std::size_t t = 0; t < m_point_sets.size(); ++ t)
       {
@@ -495,7 +495,7 @@ estimate_local_k_neighbor_scales(
   using parameters::get_parameter;
 
   typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::const_type PointMap;
-  typedef typename Point_set_processing_3::GetQueryPointMap<QueryPointRange, NamedParameters>::const_type QueryPointMap;
+  typedef typename CGAL::GetPointMap<QueryPointRange, NamedParameters, true, internal_np::query_point_t>::const_type QueryPointMap;
   typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
 
   typedef typename boost::property_traits<PointMap>::value_type Point_d;
@@ -663,7 +663,7 @@ estimate_local_range_scales(
   using parameters::get_parameter;
 
   typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::const_type PointMap;
-  typedef typename Point_set_processing_3::GetQueryPointMap<QueryPointRange, NamedParameters>::const_type QueryPointMap;
+  typedef typename CGAL::GetPointMap<QueryPointRange, NamedParameters, true, internal_np::query_point_t>::const_type QueryPointMap;
   typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
 
   typedef typename boost::property_traits<PointMap>::value_type Point_d;
diff --git a/Property_map/include/CGAL/property_map.h b/Property_map/include/CGAL/property_map.h
index 489305738c7..0dddeefd626 100644
--- a/Property_map/include/CGAL/property_map.h
+++ b/Property_map/include/CGAL/property_map.h
@@ -198,31 +198,43 @@ make_dereference_property_map(Iter)
 
 /// \ingroup PkgPropertyMapRef
 /// A `LvaluePropertyMap` property map mapping a key to itself (by reference).
+/// It is mutable if `T` is not `const` and non-mutable otherwise.
 ///
 /// \cgalModels `LvaluePropertyMap`
 template <typename T>
 struct Identity_property_map
 {
+/// \cond SKIP_IN_MANUAL
   typedef Identity_property_map<T> Self;
 
-  typedef T key_type; ///< typedef to `T`
-  typedef T value_type; ///< typedef to `T`
-  typedef const T& reference; ///< typedef to `const T&`
-  typedef boost::lvalue_property_map_tag category; ///< `boost::lvalue_property_map_tag`
+  typedef T key_type;
+  typedef T value_type;
+  typedef T& reference;
+  typedef boost::lvalue_property_map_tag category;
 
-  /// Access a property map element.
-  /// @param k a key which is returned as mapped value.
-  const value_type& operator[](const key_type& k) const { return k; }
+  value_type& operator[](key_type& k) const { return k; }
 
-  /// \name Put/get free functions
-  /// @{
-  friend reference get(const Self&, const key_type& k) { return k; }
+  friend value_type& get(const Self&, key_type& k) { return k; }
+  friend const value_type& get(const Self&, const key_type& k) { return k; }
   friend void put(const Self&, key_type& k, const value_type& v) { k = v; }
-  /// @}
+/// \endcond
 };
 
-
 /// \cond SKIP_IN_MANUAL
+template <typename T>
+struct Identity_property_map<const T>
+{
+  typedef Identity_property_map<const T> Self;
+
+  typedef T key_type;
+  typedef T value_type;
+  typedef const T& reference;
+  typedef boost::lvalue_property_map_tag category;
+
+  const value_type& operator[](key_type& k) const { return k; }
+  friend const value_type& get(const Self&, const key_type& k) { return k; }
+};
+
 template <typename T>
 struct Identity_property_map_no_lvalue
 {
@@ -399,7 +411,7 @@ struct Property_map_to_unary_function{
   {}
 
   template <class KeyType>
-  result_type
+  decltype(auto)
   operator()(const KeyType& a) const
   {
     return get(map,a);
diff --git a/Stream_support/include/CGAL/IO/Generic_writer.h b/Stream_support/include/CGAL/IO/Generic_writer.h
index 2f5d2a50d4a..b3d46b4fb77 100644
--- a/Stream_support/include/CGAL/IO/Generic_writer.h
+++ b/Stream_support/include/CGAL/IO/Generic_writer.h
@@ -54,7 +54,7 @@ public:
     m_writer.write_header(m_os, points.size(), 0, polygons.size());
     for(std::size_t i=0, end=points.size(); i<end; ++i)
     {
-      const typename boost::property_traits<PointMap>::reference p = get(point_map, points[i]);
+      const typename boost::property_traits<PointMap>::value_type& p = get(point_map, points[i]);
       m_writer.write_vertex(p.x(), p.y(), p.z());
     }
 
diff --git a/Stream_support/include/CGAL/IO/PLY/PLY_writer.h b/Stream_support/include/CGAL/IO/PLY/PLY_writer.h
index 225aa740a3f..2e8bbad3a89 100644
--- a/Stream_support/include/CGAL/IO/PLY/PLY_writer.h
+++ b/Stream_support/include/CGAL/IO/PLY/PLY_writer.h
@@ -148,7 +148,7 @@ void simple_property_write(std::ostream& stream,
                            ForwardIterator it,
                            std::pair<PropertyMap, PLY_property<std::vector<T> > > map)
 {
-  const typename PropertyMap::reference value = get(map.first, *it);
+  const typename PropertyMap::value_type& value = get(map.first, *it);
 
   if(CGAL::IO::get_mode(stream) == CGAL::IO::ASCII)
   {
diff --git a/Stream_support/include/CGAL/IO/STL.h b/Stream_support/include/CGAL/IO/STL.h
index 797376e8285..530dd5035f9 100644
--- a/Stream_support/include/CGAL/IO/STL.h
+++ b/Stream_support/include/CGAL/IO/STL.h
@@ -314,7 +314,6 @@ bool write_STL(std::ostream& os,
   PointMap point_map = choose_parameter<PointMap>(get_parameter(np, internal_np::point_map));
 
   typedef typename boost::property_traits<PointMap>::value_type             Point;
-  typedef typename boost::property_traits<PointMap>::reference              Point_ref;
   typedef typename CGAL::Kernel_traits<Point>::Kernel                       K;
   typedef typename K::Vector_3                                              Vector_3;
 
@@ -331,9 +330,9 @@ bool write_STL(std::ostream& os,
 
     for(const Triangle& face : facets)
     {
-      const Point_ref p = get(point_map, points[face[0]]);
-      const Point_ref q = get(point_map, points[face[1]]);
-      const Point_ref r = get(point_map, points[face[2]]);
+      const Point& p = get(point_map, points[face[0]]);
+      const Point& q = get(point_map, points[face[1]]);
+      const Point& r = get(point_map, points[face[2]]);
 
       const Vector_3 n = collinear(p,q,r) ? Vector_3(1,0,0) : unit_normal(p,q,r);
 
@@ -352,9 +351,9 @@ bool write_STL(std::ostream& os,
     os << "solid\n";
     for(const Triangle& face : facets)
     {
-      const Point_ref p = get(point_map, points[face[0]]);
-      const Point_ref q = get(point_map, points[face[1]]);
-      const Point_ref r = get(point_map, points[face[2]]);
+      const Point& p = get(point_map, points[face[0]]);
+      const Point& q = get(point_map, points[face[1]]);
+      const Point& r = get(point_map, points[face[2]]);
 
       const Vector_3 n = collinear(p,q,r) ? Vector_3(1,0,0) : unit_normal(p,q,r);
       os << "facet normal " << n << "\nouter loop\n";