Automated georeferencing of nighttime ISS timelapses using synthetic Earth renders, deep feature matching, and VIIRS reference data

Automated georeferencing of nighttime ISS timelapses using synthetic Earth renders, deep feature matching, and VIIRS reference data Raúl Pérez-Parrasa,b, Alejandro Sánchez de Miguela,1, Sergio Pascualb, Joaquín López-Herráizb aInstituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, Granada, 18008, Spain bUniversidad Complutense de Madrid, Faculty of Physics, Madrid, 28040, Spain Abstract Nighttime photographs acquired from the International Space Station (ISS) represent a unique resource for analysing light pollution, urban structure and nocturnal Earth dynamics. However, the lack of camera orientation metadata and the complex acquisition geometry have historically limited their scientific use. This work presents an automated pipeline for georeferencing ISS nighttime timelapses by combining physically based simulations of the ISS orbit, deep-learning feature matching between real and synthetic imagery, and refinement with VIIRS DNB satellite products. First, a global VIIRS nighttime lights texture, ISS TLE data and the camera model are used in Blender to generate synthetic views matching the time span and geometry of the real timelapse. Deep feature matching between each real frame and its corresponding synthetic render provides image-to-ground correspondences, which are projected into geographic coordinates and used as ground control points for a first georeferencing. In an optional refinement stage, georeferenced ISS frames are locally aligned to VIIRS subsets through optical flow, and the resulting displacement fields are used to correct the control points and perform a second georeferencing. The resulting workflow enables large-scale processing of ISS datasets and provides georeferenced nighttime maps suitable for ecological analysis, urban monitoring and light pollution research. Keywords: ISS imagery, VIIRS DNB, georeferencing, synthetic simulations, feature matching, optical flow, light pollution 1. Introduction Nighttime images taken from the International Space Station (ISS) represent a valuable but under-used data set to study light pollution, urban morphology, and nocturnal Earth observation (Sánchez de Miguel et al., 2019, 2021; Elvidge et al., 2017; Miller et al., 2013). Their scientific use is limited by the absence of accurate camera orientation metadata, variability in acquisition geometry, and inconsistent radiometric conditions. This work proposes an automated pipeline that combines physically based simulations of the ISS orbit, deep-learning feature matching between real and synthetic imagery, and VIIRS satellite data to georeference ISS nighttime imagery at scale. 2. Data and Methods 2.1. ISS imagery and preprocessing The input data for this work are nighttime photographs acquired from the International Space Station (ISS) and distributed through the NASA Earth Observations Laboratory (EOL) image database. For each experiment, we select a continuous sequence of images corresponding to a single ISS mission and a contiguous range of image identifiers (e.g. ISS067-E-338200 to ISS067-E-338895), which we refer to as a timelapse. The images are downloaded automatically from the EOL server using the mission name and the start and end IDs as parameters. A dedicated Python script performs the download in parallel threads, handles HTTP errors with retries, and skips files that are already present on disk so that the pipeline can be resumed without redownloading all data. Once the images have been retrieved, we extract the relevant EXIF metadata from the first and last frames of the sequence. From the first image, we read the acquisition time, the focal length of the lens and the image resolution in pixels. The acquisition time is interpreted as UTC and provides the initial time stamp for the timelapse. Together with the focal length and a predefined physical sensor size, it defines the intrinsic camera model used in the subsequent simulation and projection steps. The image width and height in pixels are later used to convert between the pixel coordinates and the angular directions in the camera frame. The last image in the sequence is used to obtain the final timestamp of the timelapse. In practice, small discrepancies may exist between the camera clock and the time used in the ISS TLE ephemerides. To account for this, we 2 introduce configurable offsets in seconds, minutes, and hours that are applied uniformly to both the start and end times. The corrected timestamps define the temporal interval over which the synthetic ISS orbit simulation is run, and ensure that the simulated frames remain synchronized with the real photographs throughout the processing chain. 2.2. Earth–ISS–camera simulation model In order to simulate ISS views that are physically consistent with the real photographs, we implement a 3D Earth–ISS–camera model in Blender (Blender Online Community, 2024), controlled from Python through the iss_simulation library. The Earth is represented as a sphere of radius RBin Blender units, which is linearly related to the physical Earth radius R⊕= 6371 km. An altitude h(in km) above the Earth surface is converted to Blender units as rB=RB1 + h R⊕,(1) so that the position of the ISS in the simulation is x=  x y z  =rB  cos φcos λ cos φsin λ sin φ  ,(2) where φand λdenote latitude and longitude in radians. The inverse mapping from Cartesian coordinates to geographic coordinates is also implemented and used when projecting pixels back onto the sphere. 2.2.1. Orbital model from TLE data The position and velocity of the ISS used in the simulation are derived from publicly available Two-Line Element (TLE) sets. For each experiment we collect all TLE files covering the acquisition period and parse them into a list of EarthSatellite objects. Given a target observation time t, we select the TLE whose epoch is closest to tand propagate the corresponding orbit using a standard SGP4 implementation (Vallado et al., 2006). This provides the geocentric position and velocity vectors of the spacecraft at time t. From the propagated state vector we compute the sub-satellite point on the Earth surface, obtaining the latitude φ(t), longitude λ(t)and altitude h(t)above the reference sphere. These quantities are used both to place the camera in the Blender scene and to define the local North–East–Up frame. 3 The velocity vector projected onto the local tangent plane is also used in the forward orientation mode to define the along-track direction. In this way, the synthetic camera poses are dynamically linked to the true ISS trajectory as described by the TLE data, ensuring temporal and geometric consistency with the real photographs. 2.2.2. Camera intrinsics and field of view The synthetic camera in Blender is defined to mimic the real ISS camera using the focal length fand a physical sensor size (ws, hs)obtained from the EXIF metadata. Assuming a pinhole model, the horizontal and vertical fields of view are θh= 2 arctan ws 2f, θv= 2 arctan hs 2f,(3) and are used to configure the Blender camera. The image resolution (W, H) in pixels is taken from the real photographs and enforced as the render resolution, so that each synthetic frame shares the same pixel grid as its corresponding ISS image. 2.2.3. Local N–E–U frame and orientation modes For every ISS position we construct a local North–East–Up (N–E–U) frame on the Earth sphere. Given the sub-satellite geographic coordinates (φ, λ)obtained from TLE propagation, the unit vectors pointing Up (U), East (E) and North (N) are defined as U= (cos φcos λ, cos φsin λ, sin φ),(4) E= (−sin λ, cos λ, 0),(5) N= (−sin φcos λ, −sin φsin λ, cos φ),(6) which form an orthonormal basis after normalisation. The nadir direction is D=−U. Two orientation conventions are considered: •North mode: the off–nadir angle (hereafter pitch) lifts the line of sight from Dtowards the geographic north N. 4 •Forward mode: the pitch is applied towards the along-track direction, obtained from the projection of the ISS velocity vector vonto the local tangent plane: F=v−(v·U)U,H0=F ∥F∥. In both conventions, the initial horizontal reference direction is H0=N (North mode) or the tangent velocity direction above (Forward mode). For a pitch angle θthe unrotated viewing direction is d0= cos θD+ sin θH0,(7) which lies on the vertical plane spanned by {D,H0}. Yaw rotation. An azimuthal rotation (yaw) by an angle ϕis applied around the Up vector U: d1=RU(ϕ)d0,u1=RU(ϕ)u0,(8) where RU(ϕ)is the Rodrigues rotation matrix RU(ϕ) = cos ϕ I +(1−cos ϕ)UUT+ sin ϕ[U]×,(9) and u0is an initial “up” direction orthogonal to d0, u0=U−(U·d0)d0,u0←u0 ∥u0∥. The rightward vector is then r1=d1×u1. Roll rotation. A roll rotation by ρis finally applied around the viewing direction d1: u= cos ρu1+ sin ρr1,(10) r=−sin ρu1+ cos ρr1.(11) 5 Camera frame mapping.. Blender defines the camera optical axis as the −Z axis. Therefore the final camera basis is zcam =−d1,xcam =r,ycam =u, and the corresponding rotation matrix is Rcam =  xcam ycam zcam  ,(12) which is converted internally into a quaternion to set the camera orientation in Blender. This formulation ensures that pitch, yaw and roll in the local N–E–U frame correspond unambiguously to physical rotations of the real ISS camera relative to the Earth’s surface and the spacecraft motion. 2.2.4. Ray casting and intersection with the Earth sphere Given the camera pose and the image resolution (W, H), each pixel (x, y) is mapped to a viewing ray in world coordinates. We first convert pixel coordinates to normalised device coordinates (NDC). Using the pixel centre convention, we define u= 2 x+ 0.5 W−1−1, v = 2 y+ 0.5 H−1−1,(13) so that u, v ∈[−1,1] by construction: the left, right, top and bottom edges correspond to u=−1,u= 1,v=−1and v= 1, respectively. This convention ensures that (0,0) corresponds to the centre of the image and that pixels are symmetrically distributed across the camera field of view. These normalised coordinates are mapped into the camera reference frame using the horizontal field of view θhand the aspect ratio aspect =W/H: dcam =  utan(θh/2) vtan(θh/2)/aspect −1  ,(14) which corresponds to the standard pinhole projection model. The resulting vector is normalised and rotated into world coordinates using the camera rotation matrix Rcam, producing the ray direction d=Rcamdcam, while the ray origin is the camera position o. 6 The intersection between the ray and the Earth sphere of radius RBis obtained by solving ∥o+td∥2=R2 B,(15) which leads to the quadratic equation at2+bt +c= 0, a =d·d, b = 2(o−c)·d, c =∥o−c∥2−R2 B,(16) with c= (0,0,0) the centre of the Earth. If the discriminant b2−4ac ≥0, the smallest positive root yields the intersection point on the sphere. This 3D point is finally converted back to geographic coordinates (latitude φ, longitude λ) and used as a ground control point for georeferencing. 2.3. Automatic camera orientation estimation The previous subsections define a forward model that, given the ISS position, camera intrinsics and a set of orientation angles, produces a synthetic nighttime view of the Earth. In practice, however, the yaw and pitch of the handheld camera are not recorded in the EXIF metadata. We therefore estimate these angles automatically by comparing a grid of simulated images against the first real frame of the timelapse. Let ψdenote the camera yaw (azimuth) and θthe off–nadir angle (pitch) in the chosen orientation convention (North or Forward mode, Section 2.2). For a given candidate pair (ψ, θ)we render a synthetic image Isim(ψ, θ)using the Earth–ISS–camera model at the acquisition time of the real image Ireal. The ISS position (latitude, longitude, altitude) and velocity are obtained from the TLE propagation at that time and kept fixed during the search. To quantify the similarity between Isim and Ireal we combine a global structural metric with a local feature-based score. First, both images are resized to a common resolution if needed and converted to grayscale, and the structural similarity index (SSIM) is computed (Wang et al., 2004): sSSIM = SSIMIsim, Ireal. Second, we run a deep feature matching model on the RGB images, based on a SuperPoint-like detector and descriptor and a learned matcher (similar to LightGlue or SuperGlue). The matcher produces a set of putative correspondences between image coordinates in Isim and Ireal, from which we extract the number of matched pairs Nm(ψ, θ). We then define a combined score S(ψ, θ) = sSSIM(ψ, θ)+α Nm(ψ, θ), 7 where αis a small weighting factor that balances the contribution of the global structural similarity and the number of deep feature matches. In practice, αis chosen so that the SSIM term dominates the score, while still favouring configurations that produce a larger number of consistent matches. The optimal yaw and pitch are obtained by maximising S(ψ, θ)over a discrete grid. We adopt a two-stage coarse-to-fine strategy. In the coarse stage, (ψ, θ)is sampled on a regular grid within user-defined ranges, ψ∈[ψmin, ψmax], θ ∈[θmin, θmax], with step sizes ∆ψcoarse and ∆θcoarse. For each grid point we render Isim(ψ, θ), compute S(ψ, θ), and retain the best-scoring pair (ψ∗, θ∗). In the fine stage, we centre a smaller search window around (ψ∗, θ∗)and repeat the evaluation with finer angular increments ∆ψfine and ∆θfine. The final estimate (ˆ ψ, ˆ θ)is taken as the maximiser of Swithin this refined window. The roll angle ρis either fixed to a small nominal value or explored over a narrow range if needed, since large rotations around the viewing axis are less common in the considered sequences and primarily affect the tilt of the horizon. For each timelapse, the estimated angles (ˆ ψ, ˆ θ, ρ)are assumed constant across all frames, reflecting the fact that the astronaut typically maintains a stable pointing during the sequence. These angles are then used in the subsequent simulation and projection steps to generate the synthetic timelapse and to derive image-to-ground control points. A simple caching mechanism stores the best angles and score for each experiment, allowing the pipeline to reuse previously estimated orientations without recomputing the search. 2.4. Synthetic timelapse simulation Once the camera intrinsics and orientation angles have been estimated, we generate a synthetic timelapse that mimics the real ISS sequence. The goal is to obtain, for each real frame, a corresponding simulated image Isim(tk) rendered from the same approximate viewpoint and with the same image resolution, using a global nighttime radiance texture derived from VIIRS DNB products. For each experiment, we first read the EXIF metadata of the first and last ISS photographs in the sequence. From these we obtain the acquisition times tstart and tend, the focal length fand the image dimensions (W, H). An optional temporal offset ∆tcan be applied to both timestamps to compensate 8 for possible camera clock drifts with respect to the TLE epoch. The resulting interval [tstart +∆t, tend +∆t]defines the time span of the synthetic timelapse. The Earth is represented in Blender as a sphere of radius RBtextured with a global nighttime lights mosaic based on VIIRS DNB data (Miller et al., 2013; Elvidge et al., 2017; Sánchez de Miguel et al., 2021). The texture is mapped in geographic coordinates as described in Section 2.2.4, so that each point on the sphere corresponds to a unique latitude–longitude pair (φ, λ). The synthetic camera is configured with the same focal length fand image size (W, H)as the real photographs, using the pinhole projection model and field-of-view formulation introduced in Section 2.2.2. The estimated yaw, pitch and roll angles (ˆ ψ, ˆ θ, ρ)from Section 2.3 are kept fixed throughout the sequence, reflecting the approximately constant pointing of the handheld camera during the timelapse. To generate the sequence, we discretise the time interval into steps of duration ∆tframe (typically between 0.5 and 1 s), obtaining a set of times {tk}such that tstart +∆t≤tk≤tend +∆t. For each time tkwe propagate the ISS orbit using the TLE-based model and compute the corresponding subsatellite point and altitude (φk, λk, hk), together with the velocity vector vk. These quantities are used to place the camera at the appropriate position above the Earth and to construct the local N–E–U frame and orientation according to the selected convention (North or Forward mode). Given (φk, λk, hk)and (ˆ ψ, ˆ θ, ρ), the simulation module renders a synthetic image Isim(tk)at resolution (W, H)using Blender’s EEVEE engine, with the background set to black and the VIIRS-based Earth texture in emission mode. The same camera intrinsics and projection model described in Section 2.2.4 are used, so that each pixel of Isim(tk)corresponds to a well-defined viewing ray intersecting the Earth sphere. The full synthetic timelapse is thus obtained as the ordered collection Isim =Isim(tk), aligned in time and geometry with the original ISS sequence. In addition to the rendered frames, we also record the ISS ground track {(φk, λk)}and optionally export a Blender scene file containing the full configuration of the timelapse. The synthetic sequence Isim is later used to establish dense correspondences with the real images and to project matched pixels onto the Earth surface in a consistent, physically-based manner. 9 VIIRS-based refinement stage described in Section 2.8, where local opticalflow corrections are estimated and a second georeferencing is optionally applied. 2.8. VIIRS-based optical-flow refinement and second georeferencing The first georeferencing step (Section 2.7) yields, for each ISS frame, a rectified image Igeo in geographic coordinates and a set of ground control points (GCPs) relating original image pixels to latitude–longitude pairs. Although this TPS-based warp already captures most of the acquisition geometry, small residual misalignments remain due to parallax, local viewing effects and the limited number of control points. To reduce these residuals, we exploit the stable VIIRS Day/Night Band (DNB) products as an external reference and estimate a dense displacement field between each georeferenced ISS image and its corresponding VIIRS subset. 2.8.1. Local ISS–VIIRS optical flow For every georeferenced ISS frame with identifier ISSxxx-E-yyyyyy, we read the RGB GeoTIFF Igeo from the geo/ directory and the matching VIIRS subset IVfrom the viirs/ directory (e.g. ISSxxx-E-yyyyyy_rect.tiff and ISSxxx-E-yyyyyy_viirs.tiff). The VIIRS subset, originally a single radiance band, is resampled to the resolution of Igeo and normalised to the range [0,255]. The ISS image is converted to grayscale and its histogram is matched to that of the VIIRS subset to reduce radiometric differences and improve the robustness of intensity-based motion estimation. Since light pollution analysis often focuses on coastal and terrestrial regions, and ocean areas contribute little structural information, we compute optical flow on a configurable sub-window of the frame. Let [x0, x1]×[y0, y1] denote this window, expressed in normalised coordinates with x, y ∈[0,1] (e.g. x∈[0,0.7],y∈[0.5,1.0] for the lower-left part of the image). The corresponding pixel bounds in the resampled images are xpix 0=⌊Wx0⌋, xpix 1=⌊Wx1⌋, ypix 0=⌊Hy0⌋, ypix 1=⌊Hy1⌋, with (W, H)the image resolution. On this cropped domain we estimate a dense optical-flow field w(x, y) = (u(x, y), v(x, y)) between the VIIRS image (taken as reference) and the histogrammatched ISS image (taken as distorted). In practice we use the Farnebäck algorithm as implemented in OpenCV (Farnebäck, 2003), with parameters 16 tuned to capture smooth, large-scale displacements while ignoring high-frequency noise. The resulting flow is stored as a two-channel array [u(x, y), v(x, y)] per pixel and saved in a NumPy file ISSxxx-E-yyyyyy_flow.npy within the flow/ directory. 2.8.2. Flow-based correction of ground control points The dense displacement field is then used to refine the ground control points generated in Section 2.6. For each image we have: •a filtered .points file containing GCPs with (mapX,mapY,sourceX,sourceY), where (sourceX,sourceY)are pixel coordinates in the original ISS frame and (mapX,mapY)are latitude and longitude from the first georeferencing; •a pixel-mapping CSV file (e.g. ISSxxx-E-yyyyyy_pixel_mapping.csv) that associates each control point with its pixel coordinates (xgeo, ygeo) in the georeferenced image; •the GeoTIFF Igeo and its geotransform, which encodes the affine mapping between pixel indices and geographic coordinates. For a given control point with rectified pixel coordinates (xgeo, ygeo), we locate the corresponding position in the cropped flow domain, compute the displacement w(xgeo, ygeo) = (u, v), and obtain a corrected pixel location x′=xgeo +u, y′=ygeo +v. If (x′, y′)falls outside the cropped region, or if no flow is available for that point, we leave the control point unchanged. Otherwise, the new pixel coordinates (x′, y′)are converted back to geographic coordinates (˜ X, ˜ Y)via the GeoTIFF geotransform, (˜ X, ˜ Y) = GeoTransform(x′, y′), yielding an updated GCP location that approximates the alignment with VIIRS. The refined control points are stored in a new set of files with suffix _corrected.points, with the same sourceX/sourceY values as the original points but updated mapX/mapY coordinates. Additional diagnostic fields record the applied pixel offsets (u, v)and a residual flag, which can be used for quality control. 17 2.8.3. Optional second georeferencing Using the corrected control points, we can optionally perform a second georeferencing pass. This step is identical in structure to the first georeferencing (Section 2.7), but replacing the original .points files by the _corrected.points versions. For each selected frame we: 1. build a new VRT with the corrected GCPs in EPSG:4326; 2. validate the GCP dispersion as in the first pass; 3. run gdalwarp with TPS interpolation to obtain an updated GeoTIFF I(2nd) geo ; 4. regenerate the pixel-mapping CSV if needed. In practice, the second pass can be applied either to the full timelapse or to a representative subset of frames (e.g. first, middle and last) to assess the benefit of the VIIRS-based correction. The refined images show improved local alignment to the VIIRS reference, particularly in regions with strong light gradients (urban cores, coastal interfaces), and provide more accurate ground control for subsequent ecological and urban analyses. 3. Results We evaluated the proposed pipeline on a set of nighttime ISS timelapses covering Europe for two different years (2017 and 2022), corresponding to Expeditions 52–53 and 67, respectively. These sequences represent some of the most challenging datasets in the ISS archive, as they include wide–angle lenses, strong geometric distortions, variable atmospheric conditions, and heterogeneous illumination. The following subsections summarise the preliminary geometric performance, the agreement with VIIRS DNB data, and the large-scale mosaics produced to date. 3.1. Geometric accuracy of the first and second georeferencing stages The first georeferencing stage, driven exclusively by physically based simulation and deep feature correspondence, typically achieves sub-kilometre consistency at regional scales. Coastlines, major roads, and large metropolitan areas appear in the correct geographic position even before refinement. After applying the VIIRS-based optical-flow correction, local distortions decrease substantially, particularly around high–contrast features such as 18 coastlines, river deltas, and compact urban cores. Although a systematic quantitative assessment is left for future work, qualitative inspection across hundreds of frames indicates: •residual planimetric errors in dense urban areas on the order of a few hundred meters; •smooth, spatially coherent displacement fields between ISS and VIIRS, consistent with the expected geometric differences between oblique and nadir views; •improved internal consistency across consecutive frames, enabling temporally stable mosaics. The main sources of residual error are low-light areas with weak spatial structure, cloud contamination, and high off-nadir angles, all of which reduce the number of reliable feature matches. 3.2. Validation using VIIRS DNB and Earth Engine mosaics To visually assess the correctness of the georeferenced frames, we compared the ISS products against VIIRS DNB composites. The ISS–VIIRS overlap reveals: •strong agreement in the spatial organisation of nighttime lights; •finer detail in ISS imagery (10–100 m/px) relative to VIIRS (500 m/px), particularly in road networks, port facilities, and intra-urban structure; •consistent colour information not available in VIIRS, enabling differentiation of lighting technologies. Using the georeferenced outputs produced by our method, calibrated mosaics were constructed in Google Earth Engine (Gorelick et al., 2017) by Sánchez de Miguel et al. (Sánchez de Miguel et al., 2019, 2021, 2025). These mosaics represent the first internally consistent, colour nighttime maps of Europe for specific years prior to the availability of colour-sensitive orbital sensors. 19 Figure 1: Nighttime mosaic of Europe for 2017 generated from ∼7,000 ISS colour photographs. The images were first georeferenced using the automated pipeline described in this work and later calibrated and mosaicked in Google Earth Engine. Black-and-white regions correspond to missing ISS coverage and are complemented with VIIRS DNB data. Different colours reflect varying lighting technologies and spectral characteristics. 3.3. Large-scale mosaics: Europe 2017 and 2022 Figure 1 shows the 2017 nighttime mosaic of Europe generated from approximately 7,000 ISS images. The mosaic leverages the georeferencing pipeline presented in this paper and was subsequently radiometrically calibrated following Sánchez de Miguel et al. (Sánchez de Miguel et al., 2021). The result preserves the large-scale distribution of nighttime lights while revealing fine-scale structure unmatched by global satellite products. Cloud-covered areas, especially over the British Isles and the Balkans, appear as regions of reduced detail. Future versions of this pipeline will incorporate automated cloud detection and removal. 3.4. Processing scalability and global potential A key outcome of this work is the full automation of the ISS georeferencing process. Once mission identifiers and image ID ranges are provided, the 20 system can process entire timelapses end-to-end without human intervention, including: 1. orientation estimation, 2. synthetic rendering, 3. deep matching, 4. GCP generation and filtering, 5. optical-flow refinement, 6. first and second georeferencing. This enables the continuous processing of large ISS archives and makes global coverage achievable. The next planned mosaics include 2012 (Expeditions 30–31) and 2024 (Expeditions 69–70), with the long-term goal of producing a global, year-by-year atlas of nighttime Earth in colour using ISS data. 4. Discussion The results obtained so far demonstrate that handheld nighttime imagery from the International Space Station can be georeferenced at scale using a combination of physically based rendering, deep feature matching and optical-flow refinement. Although the present work focuses on methodological development, the successful reconstruction of multi-thousand-image mosaics—such as the European map for 2017—indicates that the approach is robust across a wide range of acquisition geometries, resolutions and lighting conditions. A key challenge lies in the strong geometric distortions inherent to widefield ISS lenses and the lack of camera orientation metadata. The proposed synthetic-rendering strategy mitigates these limitations by embedding the camera inside a fully deterministic Earth–ISS model, allowing feature matching to operate in a physically meaningful space. The results show that even with significant optical distortion and variable exposure, the feature-based registration is stable when validated against VIIRS subsets. The optical-flow correction step further improves local alignment in heterogeneous urban regions where radiometric differences and parallax effects are more pronounced. Despite these advances, several limitations remain. First, cloud contamination is not yet handled automatically, which leads to local inconsistencies in the final mosaics, especially in northern latitudes. Second, the absence of 21 a global, high-resolution nighttime reference dataset restricts the accuracy of the optical flow, which relies on VIIRS DNB at 500 m resolution as its radiometric baseline. Third, the photometric calibration required for colourbased analysis, although available for individual frames, has not yet been integrated into the automated pipeline. Another source of uncertainty arises from the assumption of constant camera orientation across an entire timelapse. While this holds for most ISS handheld sequences, occasional small rotations can accumulate and propagate into spatial errors. Temporal smoothing or frame-to-frame orientation estimation could alleviate this. Finally, the current pipeline does not yet include global quality control metrics to automatically reject poor matches or sequences affected by severe motion blur, overexposure or partial obscuration. Nevertheless, the ability to produce continent-scale mosaics from ISS data, using only publicly available imagery and without manual annotation, highlights the robustness of the method. The approach scales linearly with the amount of imagery, opening the possibility of constructing multi-year global night-light archives from ISS photographs alone. 5. Conclusions This work introduces an automated pipeline for georeferencing nighttime ISS timelapses using synthetic Earth renders, deep-learning feature matching and VIIRS-based refinement. The method overcomes one of the main historical limitations of ISS nighttime photography: the absence of camera orientation metadata and the large geometric distortions produced by wide-field lenses. By combining deterministic simulation with data-driven matching, the pipeline produces dense ground control points without human intervention. Preliminary large-scale applications, such as the reconstruction of the 2017 and 2022 European nighttime mosaics (Sánchez de Miguel et al., 2025), confirm that the approach is operational and capable of handling tens of thousands of frames. The resulting georeferenced products maintain internal geometric consistency and are suitable for environmental and urban-lighting analyses. The method is fully reproducible, relies on open data and opensource libraries, and can run continuously across multiple timelapses, enabling the systematic processing of the full ISS nighttime archive. 22 Future developments include automatic cloud detection and removal, improved radiometric calibration for colour analyses, integration with the new Cities at Night citizen-science app, and the extension of the workflow to global coverage and to additional years (e.g. 2012 and 2024). In the long term, the combination of physically based simulation, machine learning and citizen engagement will enable the construction of a global, multi-year, highresolution atlas of nighttime lights derived exclusively from ISS observations. The work presented here demonstrates that, with appropriate modelling and automation, handheld astronaut photography can serve as a unique scientific resource for tracking the evolution of artificial light at night, complementing satellite missions and filling historical gaps in nighttime colour imaging. Acknowledgements This work received support from PLAN-B, funded by the European Union’s Horizon Europe Research and Innovation Programme (Grant Agreement No. 101135308). The author wishes to thank his thesis supervisors Alejandro Sánchez de Miguel, Sergio Pascual Ramírez and Joaquín López Herraiz for their continuous scientific guidance and support. Additional institutional support from the GUAIX group at Universidad Complutense de Madrid is gratefully acknowledged. Nighttime ISS imagery was obtained from NASA and ESA via the Earth Observations Laboratory (EOL). Funding This work was funded by PLAN-B, a project supported by the European Union’s Horizon Europe Research and Innovation Programme (Grant Agreement No. 101135308). Declaration of generative AI and AI-assisted technologies in the writing process During the preparation of this work, the author used generative AI tools to assist in structuring and editing portions of the manuscript. The author reviewed and revised all generated content and takes full responsibility for the final version of the text. 23 References Blender Online Community, 2024. Blender – a 3d modelling and rendering package. URL: https://www.blender.org. Bookstein, F.L., 1989. Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Transactions on Pattern Analysis and Machine Intelligence . DeTone, D., Malisiewicz, T., Rabinovich, A., 2018. Superpoint: Selfsupervised interest point detection and description, in: Proceedings of the CVPR Workshops. Elvidge, C.D., Baugh, K., Zhizhin, M., Hsu, F.C., 2017. Why VIIRS data are superior to DMSP for mapping nighttime lights, in: Proceedings of the Asia-Pacific Advanced Network. Farnebäck, G., 2003. Two-frame motion estimation based on polynomial expansion, in: Proceedings of the SCIA Scandinavian Conference on Image Analysis. Gorelick, N., Hancher, M., Dixon, M., et al., 2017. Google earth engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment . Luong, K., et al., 2023. Lightglue: Local feature matching at light speed. arXiv preprint arXiv:2306.02418 . Sánchez de Miguel, A., Pérez-Parras, R., Pascual, S., et al., 2025. Europe at night in real color: 2017 and 2022 mosaics from iss imagery. In preparation / internal communication. Sánchez de Miguel, A., et al., 2019. Colour remote sensing from the international space station: Mapping urban lighting. Remote Sensing of Environment . Sánchez de Miguel, A., et al., 2021. First estimation of global trends in nocturnal power emissions from VIIRS DNB data. Remote Sensing of Environment . Miller, S.D., et al., 2013. Illuminating the capabilities of the VIIRS day/night band. Remote Sensing . 24 Sarlin, P.E., DeTone, D., Malisiewicz, T., Rabinovich, A., 2020. Superglue: Learning feature matching with graph neural networks, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). Vallado, D.A., Crawford, P., Hujsak, R., Kelso, T.S., 2006. Revisiting spacetrack report #3, in: AIAA/AAS Astrodynamics Specialist Conference. Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P., 2004. Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image Processing . 25