diff --git a/e2e/features/test_get.py b/e2e/features/test_get.py
index a7d89214a..8e5b0ca37 100644
--- a/e2e/features/test_get.py
+++ b/e2e/features/test_get.py
@@ -1,7 +1,13 @@
 import pytest
-import siibra
 from e2e.util import check_duplicate
 
+import siibra
+from siibra.locations import Location, PointCloud
+from siibra.livequeries.bigbrain import WagstylProfileLoader
+from siibra.features.feature import CompoundFeature
+
+import numpy as np
+
 
 # We get all registered subclasses of Feature
 @pytest.mark.parametrize(
@@ -67,8 +73,53 @@ def test_querying_with_volume():
     # Use features with location anchors only. Because hybrid ones will also
     # employ semantic links between regions, potentially changing the result.
     region = siibra.get_region("julich 2.9", "ca1")
-    volume = region.get_regional_mask('mni152')
+    volume = region.get_regional_mask("mni152")
     profiles_region = siibra.features.get(region, "BigBrainIntensityProfile")[0]
     profiles_volume = siibra.features.get(volume, "BigBrainIntensityProfile")[0]
     # the ids will be different but the content has to be the same. Even the order.
-    assert [p.location for p in profiles_region] == [p.location for p in profiles_volume]
+    assert [p.location for p in profiles_region] == [
+        p.location for p in profiles_volume
+    ]
+
+
+WagstylProfileLoader._load()
+verts = WagstylProfileLoader._vertices
+np.random.seed(10)
+loc_queries = [
+    (
+        siibra.get_map("julich 2.9", "bigbrain")
+        .get_volume("hoc1 left")
+        .get_boundingbox(clip=True),
+        "BigBrainIntensityProfile",
+        13968,
+    ),
+    (
+        PointCloud(verts[np.random.randint(0, len(verts) - 1, 100)], space="bigbrain"),
+        "BigBrainIntensityProfile",
+        100,
+    ),
+    (
+        PointCloud(
+            np.concatenate(
+                [
+                    verts[np.random.randint(0, len(verts) - 1, 100)] + 1000,
+                    verts[np.random.randint(0, len(verts) - 1, 10)],
+                ]
+            ),
+            space="bigbrain",
+        ),
+        "BigBrainIntensityProfile",
+        10,
+    ),
+]
+
+
+@pytest.mark.parametrize("loc, feat_type, res_len", loc_queries)
+def test_query_w_locations(loc: Location, feat_type: str, res_len: int):
+    results = siibra.features.get(loc, feat_type)
+    assert len(results) > 0
+    assert results[0].data is not None
+    if len(results) == 1 and isinstance(results[0], CompoundFeature):
+        assert len(results[0]) == res_len
+    else:
+        assert len(results) == res_len
diff --git a/siibra/core/region.py b/siibra/core/region.py
index d0c5739fd..df003c6b7 100644
--- a/siibra/core/region.py
+++ b/siibra/core/region.py
@@ -124,7 +124,7 @@ def __init__(
 
     def get_related_regions(self) -> Iterable["RegionRelationAssessments"]:
         """
-        Get assements on relations of this region to others defined on EBRAINS.
+        Get assessments on relations of this region to others defined on EBRAINS.
 
         Yields
         ------
@@ -918,7 +918,7 @@ def get_peid_from_region(region: Region) -> str:
 
 def get_related_regions(region: Region) -> Iterable["RegionRelationAssessments"]:
     """
-    Get assements on relations of a region to others defined on EBRAINS.
+    Get assessments on relations of a region to others defined on EBRAINS.
 
     Parameters
     ----------
diff --git a/siibra/livequeries/bigbrain.py b/siibra/livequeries/bigbrain.py
index e5fe25cc2..a5f11b663 100644
--- a/siibra/livequeries/bigbrain.py
+++ b/siibra/livequeries/bigbrain.py
@@ -19,7 +19,7 @@
 from ..features.tabular import bigbrain_intensity_profile, layerwise_bigbrain_intensities
 from ..features import anchor as _anchor
 from ..commons import logger
-from ..locations import point, pointcloud
+from ..locations import point, pointcloud, location
 from ..core import structure
 from ..retrieval import requests, cache
 from ..retrieval.datasets import GenericDataset
@@ -123,11 +123,16 @@ def query(self, concept: structure.BrainStructure, **kwargs) -> List[bigbrain_in
         matched = concept.intersection(mesh_vertices)  # returns a reduced PointCloud with og indices as labels
         if matched is None:
             return []
+        if isinstance(matched, point.Point):
+            matched = pointcloud.from_points([matched])
         assert isinstance(matched, pointcloud.PointCloud)
+        if isinstance(concept, location.Location):
+            mesh_as_list = mesh_vertices.as_list()
+            matched.labels = [mesh_as_list.index(v.coordinate) for v in matched]
         indices = matched.labels
         assert indices is not None
         features = []
-        for i in matched.labels:
+        for i in indices:
             anchor = _anchor.AnatomicalAnchor(
                 location=point.Point(loader._vertices[i], space='bigbrain'),
                 region=str(concept),
@@ -160,12 +165,16 @@ def query(self, concept: structure.BrainStructure, **kwargs) -> List[layerwise_b
 
         loader = WagstylProfileLoader()
         mesh_vertices = pointcloud.PointCloud(loader._vertices, space='bigbrain')
-        matched = concept.intersection(mesh_vertices)  # returns a reduced PointCloud with og indices as labels
+        matched = concept.intersection(mesh_vertices)  # returns a reduced PointCloud with og indices as labels if the concept is a region
         if matched is None:
             return []
+        if isinstance(matched, point.Point):
+            matched = pointcloud.from_points([matched])
         assert isinstance(matched, pointcloud.PointCloud)
+        if isinstance(concept, location.Location):
+            mesh_as_list = mesh_vertices.as_list()
+            matched.labels = [mesh_as_list.index(v.coordinate) for v in matched]
         indices = matched.labels
-        assert indices is not None
         matched_profiles = loader._profiles[indices, :]
         boundary_depths = loader._boundary_depths[indices, :]
         # compute array of layer labels for all coefficients in profiles_left
diff --git a/siibra/locations/boundingbox.py b/siibra/locations/boundingbox.py
index c1299f1c5..019aeb0e8 100644
--- a/siibra/locations/boundingbox.py
+++ b/siibra/locations/boundingbox.py
@@ -148,11 +148,7 @@ def intersection(self, other: 'BrainStructure', dims=[0, 1, 2]):
             points_inside = [p for p in other if self.intersects(p)]
             if len(points_inside) == 0:
                 return None
-            result = pointcloud.PointCloud(
-                points_inside,
-                space=other.space,
-                sigma_mm=[p.sigma for p in points_inside]
-            )
+            result = pointcloud.from_points(points_inside)
             return result[0] if len(result) == 1 else result  # if PointCloud has single point return as a Point
 
         return other.intersection(self)
diff --git a/siibra/locations/pointcloud.py b/siibra/locations/pointcloud.py
index e2ffa9f6d..f45836087 100644
--- a/siibra/locations/pointcloud.py
+++ b/siibra/locations/pointcloud.py
@@ -120,18 +120,13 @@ def intersection(self, other: location.Location):
         if not isinstance(other, (point.Point, PointCloud, _boundingbox.BoundingBox)):
             return other.intersection(self)
 
-        intersections = [(i, p) for i, p in enumerate(self) if p.intersects(other)]
-        if len(intersections) == 0:
+        if isinstance(other, PointCloud):
+            intersecting_points = [p for p in self if p.coordinate in other.coordinates]
+        else:
+            intersecting_points = [p for p in self if p.intersects(other)]
+        if len(intersecting_points) == 0:
             return None
-        ids, points = zip(*intersections)
-        labels = None if self.labels is None else [self.labels[i] for i in ids]
-        sigma = [p.sigma for p in points]
-        intersection = PointCloud(
-            points,
-            space=self.space,
-            sigma_mm=sigma,
-            labels=labels
-        )
+        intersection = from_points(intersecting_points)
         return intersection[0] if len(intersection) == 1 else intersection
 
     @property