diff --git a/openfpm_data b/openfpm_data
index 7315141b609823336ec319370991c0a4b36e4deb..a227090455a3364dec5964b762bc2165658af189 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit 7315141b609823336ec319370991c0a4b36e4deb
+Subproject commit a227090455a3364dec5964b762bc2165658af189
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index e5416e3e84c7c5a673d78d2ebe9ec7b50c71f286..d3e359e9716e3b75080b87f777b4d5f8b6b1093f 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -3098,6 +3098,32 @@ public:
 		}
 	}
 
+	/*! \brief Remove the points
+	 *
+	 * \param box Remove all the points in the box
+	 *
+	 */
+	void removePoints(Box<dim,size_t> & box)
+	{
+		Box<dim,long int> box_i = box;
+
+		for (size_t i = 0 ; i < loc_grid.size() ; i++)
+		{
+			Box<dim,long int> bx = gdb_ext.get(i).Dbox + gdb_ext.get(i).origin;
+
+			Box<dim,long int> bout;
+			bool inte = bx.Intersect(box,bout);
+			bout -= gdb_ext.get(i).origin;
+
+			if (inte == true)
+			{
+				loc_grid.get(i).copyRemoveReset();
+				loc_grid.get(i).remove(bout);
+				loc_grid.get(i).removePoints(v_cl.getmgpuContext());
+			}
+		}
+	}
+
 	/*! \brief Move the memory from the device to host memory
 	 *
 	 */
diff --git a/src/Grid/grid_dist_id_comm.hpp b/src/Grid/grid_dist_id_comm.hpp
index d711a15bde13d901c4845698f972f6537f558565..ec045d96547d48608f483d5bb9fd266a526eeac7 100644
--- a/src/Grid/grid_dist_id_comm.hpp
+++ b/src/Grid/grid_dist_id_comm.hpp
@@ -200,7 +200,8 @@ class grid_dist_id_comm
 											  openfpm::vector<device_grid> & loc_grid,
 											  std::unordered_map<size_t,size_t> & g_id_to_external_ghost_box,
 											  const grid_sm<dim,void> & ginfo,
-											  bool use_bx_def)
+											  bool use_bx_def,
+											  size_t opt)
 	{
 		for (size_t i = 0 ; i < loc_grid.size() ; i++)
 		{loc_grid.get(i).copyRemoveReset();}
@@ -258,19 +259,24 @@ class grid_dist_id_comm
 			}
 		}
 
+		rem_copy_opt opt_ = rem_copy_opt::NONE_OPT;
+		if (opt & SKIP_LABELLING == true)
+		{opt_ = rem_copy_opt::KEEP_GEOMETRY;}
+
+
 		for (size_t i = 0 ; i < loc_grid.size() ; i++)
 		{
-			loc_grid.get(i).template removeCopyToFinalize<prp ...>(v_cl.getmgpuContext(), rem_copy_opt::PHASE1);
+			loc_grid.get(i).template removeCopyToFinalize<prp ...>(v_cl.getmgpuContext(), rem_copy_opt::PHASE1 | opt_);
 		}
 
 		for (size_t i = 0 ; i < loc_grid.size() ; i++)
 		{
-			loc_grid.get(i).template removeCopyToFinalize<prp ...>(v_cl.getmgpuContext(), rem_copy_opt::PHASE2);
+			loc_grid.get(i).template removeCopyToFinalize<prp ...>(v_cl.getmgpuContext(), rem_copy_opt::PHASE2 | opt_);
 		}
 
 		for (size_t i = 0 ; i < loc_grid.size() ; i++)
 		{
-			loc_grid.get(i).template removeCopyToFinalize<prp ...>(v_cl.getmgpuContext(), rem_copy_opt::PHASE3);
+			loc_grid.get(i).template removeCopyToFinalize<prp ...>(v_cl.getmgpuContext(), rem_copy_opt::PHASE3 | opt_);
 		}
 	}
 
@@ -1130,7 +1136,7 @@ public:
 
 		queue_recv_data_get<prp_object>(eg_box,prp_recv,prRecv_prp);
 
-		ghost_get_local<prp...>(loc_ig_box,loc_eg_box,gdb_ext,loc_grid,g_id_to_external_ghost_box,ginfo,use_bx_def);
+		ghost_get_local<prp...>(loc_ig_box,loc_eg_box,gdb_ext,loc_grid,g_id_to_external_ghost_box,ginfo,use_bx_def,opt);
 
 		for (size_t i = 0 ; i < loc_grid.size() ; i++)
 		{loc_grid.get(i).removeAddUnpackReset();}
diff --git a/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu b/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu
index 6999b460aba21c3475815ce619d16f6a8d38abbe..e0a79e8af450795f365bd148e7375091255d8e4d 100644
--- a/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu
+++ b/src/Grid/tests/sgrid_dist_id_gpu_unit_tests.cu
@@ -487,4 +487,114 @@ BOOST_AUTO_TEST_CASE( sgrid_gpu_test_conv2_test_3d )
 	BOOST_REQUIRE_EQUAL(match,true);
 }
 
+BOOST_AUTO_TEST_CASE( sgrid_gpu_test_ghost_point_remove )
+{
+	size_t sz[3] = {60,60,60};
+	periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
+
+	Ghost<3,long int> g(1);
+
+	Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	sgrid_dist_id_gpu<3,float,aggregate<float,float,float,float>> gdist(sz,domain,g,bc);
+
+	gdist.template setBackgroundValue<0>(666);
+	gdist.template setBackgroundValue<1>(666);
+	gdist.template setBackgroundValue<2>(666);
+	gdist.template setBackgroundValue<3>(666);
+
+	/////// GPU insert + flush
+
+	Box<3,size_t> box({1,1,1},{sz[0]-1,sz[1]-1,sz[2]-1});
+
+	/////// GPU Run kernel
+
+	float c = 5.0;
+
+	typedef typename GetAddBlockType<decltype(gdist)>::type InsertBlockT;
+
+	gdist.addPoints(box.getKP1(),box.getKP2(),
+			        [] __device__ (int i, int j, int k)
+			        {
+						return true;
+			        },
+			        [c] __device__ (InsertBlockT & data, int i, int j, int k)
+			        {
+			        	data.template get<0>() = c + i + j + k;
+			        	data.template get<1>() = c + 1000 + i + j + k;
+			        }
+			        );
+
+	gdist.template flush<smax_<0>,smax_<1>>(flush_type::FLUSH_ON_DEVICE);
+
+	gdist.template ghost_get<0,1>(RUN_ON_DEVICE);
+
+	// Remove the right side of the points
+	Box<3,size_t> bxR({59,0,0},{59,59,59});
+	gdist.removePoints(bxR);
+
+	// Remove the right side of the points
+	Box<3,size_t> bxT({0,0,59},{59,59,59});
+	gdist.removePoints(bxT);
+
+	// Remove the right side of the points
+	Box<3,size_t> bxD({0,59,0},{59,59,59});
+	gdist.removePoints(bxD);
+
+	gdist.template ghost_get<0,1>(RUN_ON_DEVICE);
+
+	gdist.deviceToHost<0>();
+	gdist.write_debug("Test_out");
+
+/*	for (int i = 0 ; i < 10 ; i++)
+	{
+		gdist.template ghost_get<0,1>(RUN_ON_DEVICE);
+	}
+
+	// Now run the convolution
+
+	typedef typename GetCpBlockType<decltype(gdist),0,1>::type CpBlockType;
+
+	gdist.template conv2<0,1,2,3,1>({2,2,2},{(int)sz[0]-2,(int)sz[1]-2,(int)sz[2]-2},[] __device__ (float & u_out, float & v_out, CpBlockType & u, CpBlockType & v,int i, int j, int k){
+		u_out = u(i+1,j,k) - u(i-1,j,k) + u(i,j+1,k) - u(i,j-1,k) + u(i,j,k+1) - u(i,j,k-1);
+		v_out = v(i+1,j,k) - v(i-1,j,k) + v(i,j+1,k) - v(i,j-1,k) + v(i,j,k+1) - v(i,j,k-1);
+	});
+
+	gdist.deviceToHost<0,1,2,3>();
+
+	// Now we check that ghost is correct
+
+	auto it3 = gdist.getSubDomainIterator({2,2,2},{(int)sz[0]-2,(int)sz[1]-2,(int)sz[2]-2});
+
+	bool match = true;
+
+	while (it3.isNext())
+	{
+		auto p = it3.get();
+
+		auto p_xp1 = p.move(0,1);
+		auto p_xm1 = p.move(0,-1);
+		auto p_yp1 = p.move(1,1);
+		auto p_ym1 = p.move(1,-1);
+		auto p_zp1 = p.move(2,1);
+		auto p_zm1 = p.move(2,-1);
+
+		float sub1 = gdist.template get<2>(p);
+		float sub2 = gdist.template get<3>(p);
+
+		if (sub1 != 6.0 || sub2 != 6.0)
+		{
+			std::cout << sub1 << "  " << sub2 << std::endl;
+			std::cout << gdist.template get<0>(p_xp1) << "   " << gdist.template get<0>(p_xm1) << std::endl;
+			std::cout << gdist.template get<1>(p_xp1) << "   " << gdist.template get<1>(p_xm1) << std::endl;
+			match = false;
+			break;
+		}
+
+		++it3;
+	}
+
+	BOOST_REQUIRE_EQUAL(match,true);*/
+}
+
 BOOST_AUTO_TEST_SUITE_END()