Run a command at every step of a relaxed scan
Guess what: It is possible (in your specific case).

I was wrong. At least a little bit. Most of what is described below is true for the general purpose, and I will not change that part of the answer.

However, in the specific case where you only want to compute the volume, you are able to do that because I think it is part of the population analysis module. You can actually request a population analysis at every step of the optimisation, but be careful what you wish for (running NBO will take forever). That way the optimisation will be a little slower, but that is the most convenient way. Therefore the solution to your very specific problem is adding the following line to the route section:


pop=always volume


Be aware, that the output will be an even bigger nightmare than before.

Here is a complete input of the $\ce{H3B-NH3}$ I have used previously.


%chk=df-bp86-d3.def2svp.chk
%nproc=4
%mem=16000MB
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
pop=always                         ! Perform a population analysis at every step
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP
Scan

0 1
  N       0.00000        0.00000       -0.80852
  H       0.92634        0.00000       -1.36025
  H      -0.46317       -0.80223       -1.36025
  H      -0.46317        0.80223       -1.36025
  B       0.00000        0.00000        0.19148
  H       0.52348        0.90670        0.85928
  H      -1.04696        0.00000        0.85928
  H       0.52348       -0.90670        0.85928

B 1 5 S 5 0.2



For comparison (I'm not 100% sure I used the right values):


Energies
Step1 SCF Done:  E(RB-P86) =  -82.7531561505     A.U. after    8 cycles
Step2 SCF Done:  E(RB-P86) =  -83.0511089272     A.U. after    7 cycles
Step3 SCF Done:  E(RB-P86) =  -83.1459883448     A.U. after    7 cycles
Step4 SCF Done:  E(RB-P86) =  -83.1661412429     A.U. after    8 cycles
Step5 SCF Done:  E(RB-P86) =  -83.1606918387     A.U. after    8 cycles
Step6 SCF Done:  E(RB-P86) =  -83.1485268431     A.U. after    8 cycles

Volumes
Step1 Molar volume =  491.032 bohr**3/mol ( 43.819 cm**3/mol)
Step2 Molar volume =  415.148 bohr**3/mol ( 37.047 cm**3/mol)
Step3 Molar volume =  401.060 bohr**3/mol ( 35.790 cm**3/mol)
Step4 Molar volume =  406.785 bohr**3/mol ( 36.301 cm**3/mol)
Step5 Molar volume =  514.428 bohr**3/mol ( 45.907 cm**3/mol)
Step6 Molar volume =  453.278 bohr**3/mol ( 40.450 cm**3/mol)


The energies differ very slightly only (from the below approach), but the volumes do differ at some points. Now that could be, that I just used the wrong population analysis, or something else. I don't have the time to investigate. (Leave a comment or edit if you find out.)



I've again looked into that, but without some serious hacking via the external keyword interface, I don't think it is possible. I also couldn't find a suitable iOP to control what gets punched to checkpoint file and when, and then use the result from that for a property/analysis run.
The punch keyword is equally disappointing, only outputting the last state.

While the manual suggests that you can access a given step of a scan in the checkpoint file, I found this to be not true. The documentation is fuzzy on that part and it also warns, that not everything is stored. Since it is a binary file, there is also no real way of checking what is contained.
I guess the real reason for the existence of that option is to restart a calculation that went rogue with different parameters from an intermediate calculation.

I was equally unable to find a way to set a given redundant coordinate to a specific value and then freezing it in order to perform a partial optimisation at that point, and then manipulate that coordinate while reading in the guess and perform another calculation.
Even with the new and generally awesome way to define general internal coordinates in Gaussian 16 (GIC) I was unsuccessful. 

From all that trying, I came up with the probably easiest (=lazy) way (not computationally most efficient) to do what you want is to perform the scan, extract the (p'optimised) coordinates and run a series of single point calculations.  

A slightly less (probably) computationally heavy approach would be to perform an initial scan on a very low level (like pm6, if your molecule is well-behaved), and the do p'optimisations on the resulting steps. Usually pm6 is able to offer reasonable starting geometries, which in many other cases may also help you reduce computational effort. So this is probably easy and quite efficient.
There are various options to extract geometries from a scan, the following question deals exactly with that: Extract all structures of Gaussian 09 molecular dynamics calculation using babel?

There certainly are more efficient ways to do this kind of task, but in the end, most of them require manual labour.



N.B.: I cannot believe it, but the pm6/popt approach is actually a halfway decent application of compound jobs. Although I have to admit I have not tried it yet.
It works as the following example demonstrates. (I used Gaussian 16 Rev. A.03, but it should work with Gaussian 09 Rev D.01, too. Maybe even earlier.)

I set up an initial scan on the pm6 level of theory for $\ce{H3B-NH3}$ and scanned the $\ce{B-N}$ bond:


%chk=pm6.scan.chk
%nproc=2
%mem=8000MB
#P PM6
OPT(MaxCycle=100)
SYMMETRY(loose)
GEOM(ModRedundant)
Volume

title

0 1
 N     0.000000     0.000000    -0.764994
 H     1.019690     0.000000    -1.164250
 H    -0.509845    -0.883078    -1.164250
 H    -0.509845     0.883078    -1.164250
 B    -0.000000     0.000000     0.235006
 H     0.509845     0.883078     0.634262
 H    -1.019691     0.000000     0.634262
 H     0.509845    -0.883078     0.634262

B 1 5 S 5 0.2



I have extracted the pre-optimised geometries with Chemcraft and added a header section and made use of the --Link1-- separator. After the first optimisation I read the MO for a slight speedup. I'm sure some of this can be further automatised, but for the sake of demonstrating, I did it by hand. The commented input:


%chk=df-bp86-d3.def2svp.rescan.chk
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP 
Step 1

0 1
  N       0.00000        0.00000       -0.80852
  H       0.92634        0.00000       -1.36025
  H      -0.46317       -0.80223       -1.36025
  H      -0.46317        0.80223       -1.36025
  B       0.00000        0.00000        0.19148
  H       0.52348        0.90670        0.85928
  H      -1.04696        0.00000        0.85928
  H       0.52348       -0.90670        0.85928

B 1 5 F

--Link1--
%chk=df-bp86-d3.def2svp.rescan.chk
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
guess(read)                        ! Use the MO from the previous step
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP
Step 2

0 1
  N       0.00000        0.00000       -0.90852
  H       0.94043        0.00000       -1.36414
  H      -0.47021       -0.81444       -1.36414
  H      -0.47021        0.81444       -1.36414
  B       0.00000        0.00000        0.29148
  H       0.53853        0.93276        0.86317
  H      -1.07706        0.00000        0.86317
  H       0.53853       -0.93276        0.86317

B 1 5 F

--Link1--
%chk=df-bp86-d3.def2svp.rescan.chk
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
guess(read)                        ! Use the MO from the previous step
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP
Step 3

0 1
  N       0.00000        0.00000       -0.99853
  H       0.95023        0.00000       -1.38525
  H      -0.47511       -0.82292       -1.38525
  H      -0.47512        0.82292       -1.38525
  B       0.00000        0.00000        0.40147
  H       0.55372        0.95907        0.87762
  H      -1.10743        0.00000        0.87762
  H       0.55372       -0.95907        0.87762

B 1 5 F

--Link1--
%chk=df-bp86-d3.def2svp.rescan.chk
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
guess(read)                        ! Use the MO from the previous step
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP
Step 4

0 1
  N       0.00000        0.00000       -1.07752
  H       0.95590        0.00000       -1.42049
  H      -0.47795       -0.82783       -1.42049
  H      -0.47795        0.82783       -1.42049
  B       0.00000        0.00000        0.52248
  H       0.56746        0.98288        0.89885
  H      -1.13493        0.00000        0.89885
  H       0.56746       -0.98288        0.89885

B 1 5 F

--Link1--
%chk=df-bp86-d3.def2svp.rescan.chk
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
guess(read)                        ! Use the MO from the previous step
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP
Step 5

0 1
  N       0.00000        0.00000       -1.14730
  H       0.95649        0.00000       -1.46941
  H      -0.47824       -0.82834       -1.46941
  H      -0.47824        0.82834       -1.46941
  B       0.00000        0.00000        0.65270
  H       0.57849        1.00197        0.92763
  H      -1.15698        0.00000        0.92763
  H       0.57849       -1.00197        0.92763

B 1 5 F

--Link1--
%chk=df-bp86-d3.def2svp.rescan.chk
#P BP86/def2SVP/W06                ! Density Functional Theory Calculation
DenFit                             ! Use density fitting
empiricaldispersion=GD3            ! Use Grimme Dispersion
guess(read)                        ! Use the MO from the previous step
opt(MaxCycle=100,Loose,ModRed)     ! Use more optimisation cycles (Loose only for speed)
scf(xqc,MaxConventionalCycle=500)  ! If necessary, resort to quadratic convergence
int(ultrafinegrid)                 ! Larger Grid
scrf(pcm,solvent=water)            ! Use solvent
gfinput gfoldprint iop(6/7=3)      ! For molden
symmetry(loose)                    ! Loosen symmetry requirements
Volume                             ! Report vdW volume

Water H3B-NH3 DF-BP86-D3(PCM)/def2-SVP
Step 6

0 1
  N       0.00000        0.00000       -1.21413
  H       0.95402        0.00000       -1.53258
  H      -0.47701       -0.82621       -1.53258
  H      -0.47701        0.82621       -1.53258
  B       0.00000        0.00000        0.78587
  H       0.58528        1.01374        0.96868
  H      -1.17057        0.00000        0.96868
  H       0.58528       -1.01374        0.96868

B 1 5 F



For comparison:
Energies


Step1 SCF Done:  E(RB-P86) =  -82.7531561507     A.U. after    8 cycles
Step2 SCF Done:  E(RB-P86) =  -83.0511087360     A.U. after    8 cycles
Step3 SCF Done:  E(RB-P86) =  -83.1459860378     A.U. after    7 cycles
Step4 SCF Done:  E(RB-P86) =  -83.1661471840     A.U. after    8 cycles
Step5 SCF Done:  E(RB-P86) =  -83.1606709617     A.U. after    8 cycles
Step5 SCF Done:  E(RB-P86) =  -83.1484921220     A.U. after    8 cycles


Volumes


Step1 Molar volume =  459.690 bohr**3/mol ( 41.022 cm**3/mol)
Step2 Molar volume =  479.017 bohr**3/mol ( 42.747 cm**3/mol)
Step3 Molar volume =  401.060 bohr**3/mol ( 35.790 cm**3/mol)
Step4 Molar volume =  406.785 bohr**3/mol ( 36.301 cm**3/mol)
Step5 Molar volume =  570.344 bohr**3/mol ( 50.897 cm**3/mol)
Step6 Molar volume =  498.606 bohr**3/mol ( 44.495 cm**3/mol)


Animation