Geodetically constrained interseismic interplate coupling has been widely used to assess seismic potential in subduction zones. Modeling interseismic deformation is challenging, as it involves interplate coupling and often ignores continental internal deformation processes. We present a novel methodology to jointly estimate interplate coupling along with upper plate rigid motion and surface strain, constrained by GNSS‐derived velocities. We use a least squares inversion with a spatially variable Equal Posterior Information Condition Tikhonov regularization, accounting for observational and elastic structure uncertainties. Our modeling reveals three megathrust regions with high tsunamigenic earthquake potential located within the Atacama Seismic Gap (Chile). This study indicates the presence of a downdip segmentation located just above the 1995 (Mw8.0) Antofagasta earthquake rupture, raising concerns for the potential of tsunamigenic earthquake occurrence at shallower depths. Additionally, we show that surface motion is dominated by strain, with rather negligible rigid motion, dismissing the rigid Andean microplate model typically assumed in previous studies. Plain Language Summary We present a novel methodology to investigate the degree of coupling (a measure of the earthquake potential) along the plate interface in subduction zones. Here, we infer deformation and rotation of the continental plate together with the degree of coupling, all using Global Positioning System (GPS) measurements of continental surface deformation. We apply this approach for the subduction margin formed by the convergence of Nazca and South‐American plates, and located between the Chilean cities of Antofagasta and La Serena. We find three regions with high seismic potential, raising concerns for future occurrence of large tsunamigenic earthquakes. The first region is located westward of the rupture area of the 1995 (Mw8.0) Antofagasta earthquake, the second region extends for 200 km between the cities of Taltal and Copiapo and the third region extends for 100 km between the cities of Vallenar and La Serena. Our results highlight the importance of estimating continental strain jointly with coupling models. Key Points Improved subduction strain accumulation model between 23° and 30°S taking advantage of newly developed inversion schemes Negligible rigid Andean block motion and fore‐arc extension dominates the Atacama Region (23°–30°S) Above the 1995 Mw8.0 earthquake, the megathrust remains coupled with high seismic potential