  Nonlinear electromagnetic response functions have reemerged as a crucial tool
for studying quantum materials. Most attention has been paid to responses to
spatially uniform electric fields, relevant to optical experiments in
conventional materials. However, magnetic and magnetoelectric phenomena are
naturally connected by responses to spatially varying electric fields due to
Maxwell's equations. Furthermore, in moir\'{e} materials, characteristic
lattice scales are much longer, allowing spatial variation of optical electric
fields to potentially have a measurable effect in experiments. To address these
issues, we develop a formalism for computing linear and nonlinear responses to
spatially inhomogeneous electromagnetic fields. Starting with the continuity
equation, we derive an expression for the second-quantized current operator
that is manifestly conserved and model independent. Crucially, our formalism
makes no assumptions on the form of the microscopic Hamiltonian and so is
applicable to model Hamiltonians derived from tight-binding or ab initio
calculations. We then develop a diagrammatic Kubo formalism for computing the
wavevector dependence of linear and nonlinear conductivities, using Ward
identities to fix the value of the diamagnetic current order-by-order in the
vector potential. We apply our formula to compute the magnitude of the Kerr
effect at oblique incidence for a model of a moir\'{e} Chern insulator and
demonstrate the experimental importance of spatially inhomogeneous fields in
these systems. We further show how our formalism allows us to compute the
(orbital) magnetic multipole moments and magnetic susceptibilities in
insulators. We next use our formalism to compute the second-order transverse
response to spatially varying transverse electric fields in our moir\'{e} Chern
insulator model, with an eye towards the next generation of experiments in
these systems.
