diff --git a/xcp_d/interfaces/workbench.py b/xcp_d/interfaces/workbench.py
index 90de8b247..f809db9c9 100644
--- a/xcp_d/interfaces/workbench.py
+++ b/xcp_d/interfaces/workbench.py
@@ -1168,3 +1168,71 @@ class CiftiCreateDenseFromTemplate(WBCommand):
     input_spec = _CiftiCreateDenseFromTemplateInputSpec
     output_spec = _CiftiCreateDenseFromTemplateOutputSpec
     _cmd = "wb_command -cifti-create-dense-from-template"
+
+
+class _CiftiChangeMappingInputSpec(CommandLineInputSpec):
+    """Input specification for the CiftiCreateDenseFromTemplate command."""
+
+    data_cifti = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=0,
+        desc="The cifti file to use the data from.",
+    )
+    direction = traits.Enum(
+        "ROW",
+        "COLUMN",
+        mandatory=True,
+        argstr="%s",
+        position=1,
+        desc="Which mapping to parcellate (integer, ROW, or COLUMN)",
+    )
+    cifti_out = File(
+        name_source=["label"],
+        name_template="converted_%s.dscalar.nii",
+        keep_extension=False,
+        argstr="%s",
+        position=2,
+        desc="The output cifti file.",
+    )
+    scalar = traits.Bool(
+        False,
+        usedefault=True,
+        argstr="-scalar",
+        position=3,
+        desc="Set the mapping to scalar",
+    )
+
+
+class _CiftiChangeMappingOutputSpec(TraitedSpec):
+    """Output specification for the CiftiCreateDenseFromTemplate command."""
+
+    cifti_out = File(exists=True, desc="output CIFTI file")
+
+
+class CiftiChangeMapping(WBCommand):
+    """Convert to scalar, copy mapping, etc.
+
+    Take an existing cifti file and change one of the mappings.
+    Exactly one of -series, -scalar, or -from-cifti must be specified.
+    The direction can be either an integer starting from 1, or the strings 'ROW' or 'COLUMN'.
+
+    Examples
+    --------
+    >>> ccdft = CiftiChangeMapping()
+    >>> ccdft.inputs.data_cifti = "tpl-fsLR_atlas-Gordon_den-32k_dseg.dlabel.nii"
+    >>> ccdft.inputs.direction = "ROW"
+    >>> ccdft.inputs.cifti_out = "out.dscalar.nii"
+    >>> ccdft.inputs.scalar = True
+    >>> ccdft.cmdline
+    wb_command -cifti-change-mapping \
+        tpl-fsLR_atlas-Gordon_den-32k_dseg.dlabel.nii \
+        ROW \
+        out.dscalar.nii \
+        -scalar
+    """
+
+    input_spec = _CiftiChangeMappingInputSpec
+    output_spec = _CiftiChangeMappingOutputSpec
+    _cmd = "wb_command -cifti-change-mapping"
diff --git a/xcp_d/workflows/connectivity.py b/xcp_d/workflows/connectivity.py
index c295cfcbc..cd3d16a19 100644
--- a/xcp_d/workflows/connectivity.py
+++ b/xcp_d/workflows/connectivity.py
@@ -17,6 +17,7 @@
 )
 from xcp_d.interfaces.nilearn import IndexImage
 from xcp_d.interfaces.workbench import (
+    CiftiChangeMapping,
     CiftiCreateDenseFromTemplate,
     CiftiParcellateWorkbench,
 )
@@ -170,6 +171,8 @@ def init_load_atlases_wf(
         # fmt:on
 
     else:
+        # Add empty vertices to atlas for locations in data, but not in atlas
+        # (e.g., subcortical regions for cortex-only atlases)
         resample_atlas_to_data = pe.MapNode(
             CiftiCreateDenseFromTemplate(),
             name="resample_atlas_to_data",
@@ -185,6 +188,26 @@ def init_load_atlases_wf(
         ])
         # fmt:on
 
+        # Change the atlas to a scalar file.
+        convert_to_dscalar = pe.MapNode(
+            CiftiChangeMapping(
+                direction="ROW",
+                scalar=True,
+                cifti_out="atlas.dscalar.nii",
+            ),
+            name="convert_to_dscalar",
+            mem_gb=mem_gb,
+            n_procs=omp_nthreads,
+            iterfield=["data_cifti"],
+        )
+        # fmt:off
+        workflow.connect([
+            (resample_atlas_to_data, convert_to_dscalar, [("cifti_out", "data_cifti")]),
+        ])
+        # fmt:on
+
+        # Convert atlas from dlabel to pscalar format.
+        # The pscalar version of the atlas is later used for its ParcelAxis.
         parcellate_atlas = pe.MapNode(
             CiftiParcellateWorkbench(
                 direction="COLUMN",
@@ -200,7 +223,7 @@ def init_load_atlases_wf(
         # fmt:off
         workflow.connect([
             (atlas_file_grabber, parcellate_atlas, [("atlas_file", "atlas_label")]),
-            (resample_atlas_to_data, parcellate_atlas, [("cifti_out", "in_file")]),
+            (convert_to_dscalar, parcellate_atlas, [("cifti_out", "in_file")]),
             (parcellate_atlas, outputnode, [("out_file", "parcellated_atlas_files")]),
         ])
         # fmt:on