## 3 Materials and Methods
### O-FDTD Simulation of LoRa Coverage  
The cityscape LoRa propagation is simulated using the O-FDTD method [13] via the (c) 2020 WaveBox toolbox for MATLAB, Supplementary Materials S1. O-FDTD describes the propagation of electric field waves by modifying the Lorentz Oscillator Model (LOM). Similarly to the LOM, which considers the oscillation of atom-bound electrons driven by an incident electric field, O-FDTD uses a mesh of coupled oscillators whose orthogonal displacement relates to the electric field amplitude. However, O-FDTD's oscillators are not electrons but abstract electric field carriers evenly distributed across space, including vacuum. Each oscillator is connected to its mesh neighbors by a refractive index-dependent coupling strength that defines the electric field transfer rate and, consequently, the wave propagation speed. The electric field \(E\) time propagation is obtained using the leap-frog backward differentiation method  
\[E_{i}=E_{i-1}+\frac{\text{d}E_{i-1}}{\text{d}t}\delta t+\frac{1}{2}\frac{\text {d}^{2}E_{i-1}}{\text{d}t^{2}}\delta t^{2}, \tag{3}\]  
where \(t\) is the time, \(\delta t\) is the temporal resolution, and the index \(i\) refers to the current time step. The single motion equation is formulated as:  
\[\frac{\text{d}^{2}E_{i-1}}{\text{d}t^{2}}=\frac{8}{3n^{2}\delta t^{2}}\Big{(}E _{i-1}^{neigh}-E_{i-1}\Big{)}-\frac{\sqrt{2}c}{n}\frac{4\pi\kappa}{\lambda} \frac{\text{d}E_{i-1}}{\text{d}t}, \tag{4}\]  
where \(n\) and \(\kappa\) are the real and imaginary parts of the refractive index, \(E^{neigh}\) is the neighborhood electric field, and \(c\) and \(\lambda\) are the vacuum speed of light and wavelength. The first term in Equation (4) is analogous to the restoring spring motion as a function of the neighborhood electric field and accounts for the wave propagation and interference effects. The second term of Equation (4) relates to the energy damping or absorption. Using this simple formulation, O-FDTD allows simulating 2D maps of material-dependent electric field wave speed and absorption. The refractive index contrast induces signal reflections at the building walls.  
The (c) 2020 WaveBox toolbox provides simple design functions and allows importing Google maps [51] data directly from a bitmap image. First, the map is converted into a vector format using Inkscape. Using Inkscape's brightness-threshold bitmap-tracing function and minor manual editing, the city map is cleaned from features, such as text, icons, and roads, so that only building walls remain. Building walls are (arbitrarily) set as blue and empty space as black. We considered a 600 \(\times\) 600 m\({}^{2}\) map surrounded by 0.60 m margins used for the perfectly matched layers (PML) boundary conditions required by the O-FDTD algorithm and exported the final map as a bitmap *.png file. Air and the stone walls were modeled in WaveBox using refractive indices of \(n_{\text{air}}\) and \(n_{\text{wall}}\), respectively (deduced in Section 4.1). The simulation parameters are summarized in Table 2. For the sake of simplicity, the simulation can be carried out in the reverse of the experimental configuration, i.e., placing a dipole source at the experimental position of the LoRa gateway (receiver), located at the center of the simulation space. The O-FDTD simulation yields a 2D map of the signal power, which we converted into RSSI values given in dBm, representing the signal propagation across the cityscape. The coverage area is defined using the calculated (Section 4.1) receiver sensitivity as a threshold in the false color map. Detectable and undetectable values are shown in color and grayscales, respectively. The simulation of a 600 \(\times\) 600 m\({}^{2}\) area took about 3 h to complete. To visualize the simulation results described in Section 4.1, we overlayed a transparent vector-image of the city map onto the obtained RSSI values using Inkscape.