When modeling global satellite data to recover a planetary magnetic or gravitational potential field and evaluate it elsewhere, the method of choice remains their analysis in terms of spherical harmonics. When only regional data are available, or when data quality varies strongly with geographic location, the inversion problem becomes severely ill-posed. In those cases, adopting explicitly local methods is to be preferred over adapting global ones (e.g., by regularization). Here, we develop the theory behind a procedure to invert for planetary potential fields from vector observations collected within a spatially bounded region at varying satellite altitude. Our method relies on the construction of spatiospectrally localized bases of functions that mitigate the noise amplification caused by downward continuation (from the satellite altitude to the planetary surface) while balancing the conflicting demands for spatial concentration and spectral limitation. The 'altitude-cognizant' gradient vector Slepian functions (AC-GVSF) were first employed in a preceding paper. They enjoy a noise tolerance under downward continuation that is much improved relative to the 'classical' gradient vector Slepian functions (CL-GVSF), which do not factor satellite altitude into their construction. Furthermore, venturing beyond the realm of their first application, in the present article we extend the theory to being able to handle both internal and external potential-field estimation. Solving simultaneously for internal and external fields in the same setting of regional data availability reduces internal-field artifacts introduced by downward-continuing unmodeled external fields, as we show with numerical examples. We explain our solution strategies on the basis of analytic expressions for the behavior of the estimation bias and variance of models for which signal and noise are uncorrelated, (essentially) spaceand bandlimited, and spectrally (almost) white. The AC-GVSF are optimal linear combinations of vector spherical harmonics. Their construction is not altogether very computationally demanding when the concentration domains (the regions of spatial concentration) have circular symmetry, e.g., on spherical caps or rings -even when the spherical-harmonic bandwidth is large. Data inversion proceeds by solving for the expansion coefficients of truncated function sequences, by least-squares analysis in a reduced-dimensional space. Hence, our method brings high-resolution regional potential-field modeling from incomplete and noisy vector-valued satellite data within reach of contemporary desktop machines.