PARKING OCCUPANCY ESTIMATION ON SENTINEL-1 IMAGES

This paper presents a method to estimate the occupancy ratio of parkings from SAR satellite images. The algorithm takes as input a series of Sentinel-1 images along with a mask indicating where the parking is located and returns for each image an occupancy ratio. The method is generic and can easily be extended. We validate our results in two parts. First, we have created a dataset of Sentinel-1 GRD image time series where each image is associated to a ground truth parking occupancy ratio. This ground truth is estimated thanks to a surveillance camera that permanently films and records the parking. We observe a strong correlation between the estimated occupancy rate and the ground truth occupancy rate. Secondly, we estimate the occupancy ratio of the 250 largest retail parkings in France from January 2018 to April 2020. We observe that weekly and seasonal patterns are consistent with consumer and economic trends. Parking occupancy estimations also plummet during the COVID-19 containment measures.


INTRODUCTION
Automatic parking occupancy monitoring can greatly help public and private institutions make better informed decisions. Monitoring hospitals parkings can be used as an early indicator of disease outbreak as described in (Butler et al., 2014) and (Nsoesie et al., 2015). It can help cities determine where to allocate parking space and how to define parkings pricing schemes (Cats et al., 2016). Monitoring stores parkings can be a proxy for estimating their activity and revenue (Partnoy, 2019).
Thus, estimating parking occupancy has been an important topic of research these last years. (Grodi et al., 2016) use sensors installed on each parking spot. (Aryandoust et al., 2019) use Uber Movement travel time data to infer parking occupancy. Although these methods are accurate, they require the installation of measuring devices in parkings or cars. They can therefore be costly to use at large scale, and require that either drivers or parkings share their data.
Others use therefore cameras or remote sensing devices. (Amato et al., 2017) and (Tȃtulea et al., 2019) use surveillance cameras and computer vision techniques to classify the occupancy of each parking spot. (Paidi, Fleyeh, 2019) use thermal cameras to detect cars. Vehicle detectors, both for standard images (Behrendt, 2019) or aerial / satellite images (Xia, et al., 2018, Drouyer, de Franchis, 2019, Drouyer, 2020b, can be used to estimate parking occupancy. Parking occupancy estimation can also be done on lower resolution optical satellites images by measuring color and gradient magnitude to differentiate occupied and unoccupied areas (Drouyer, 2020a).
Compared to other methods, using satellite images offers the advantage of not having to install any measuring device or camera near parkings. Moreover, recent satellite constellations such as Sentinel-1 allow to get free images of any area with sub- * Corresponding author: sebastien@drouyer.com weekly revisit. A significant advantage of SAR satellites such as Sentinel-1 compared to optical satellites is that the detection can be done even when there are clouds or adverse atmospheric or meteorological conditions. We present in this paper a method that estimates the occupancy ratio of open car parks from time series of Sentinel-1 images. The method takes as input a series of Sentinel-1 images along with a mask indicating where the parking is located. It returns for each image an occupancy ratio comprised between 0% and 100%.
To validate our approach, we have created a dataset of image series. We were able to get a reasonably accurate ground truth occupancy rate for one parking site (that we will call validation site) as the parking was constantly filmed by a surveillance camera and we had access to records on several months. We also have image series of 250 of the largest retail parkings in France. We don't have ground truth associated to them, but we observe that our obtained weekly and seasonal patterns are consistent with consumer and economic trends.

OBSERVATIONS
Several images of the validation site are displayed in Figure 3. We suppose in this section that we have first handled common issues that appear when dealing with SAR satellite images as described in Section 3: all images have been calibrated and orthorectified.
A general observation that can be made is that areas filled by cars tend to be much brighter. This can be explained both because of the geometry of cars and because metallic objects tend to be more reflective. This can be observed on both VV and VH polarisation channels of the SAR image, although in different intensity levels and with slight spatial differences. A notable difference is the occasional presence of sinc reflections in the VV channel, as shown in Figure 1. This makes the VV channel less easily exploitable than the VH channel.
The Sentinel-1 satellites acquire images from both the ascending and descending passes of each orbit. Images taken from descending passes are acquired early in the morning (around 6 am local time in the overflown area), and those taken from ascending passes are acquired late in the afternoon (around 6 pm). As retail parkings are empty early in the morning, only images taken from ascending passes are relevant. Images taken from ascending and descending passes tend to feature the same area in different ways, so a separated processing might be needed if one wants to use both passes. Even if restricted to a single pass, some areas can be acquired with two different incidence angles, making the image aspect vary over time. However this variation is more limited than the variation between ascending and descending passes, as shown in Figure 2. . Temporal median images of an empty parking. Images taken from ascending and descending passes were sorted apart, as well as images taken with two different incidence angles for the ascending pass.
Parkings regularly contain reflective objects / materials other than cars. These reflective objects tend to remain fixed: trees, buildings, or metallic structures. It is therefore likely that some white spots always show up at a constant position on the parking. Figure 2 shows the effect of buildings (upper left) and vegetation (lower right).
Finally, parkings don't fill in a random pattern: cars tend to park near important structures first (nearest to shops or restaurants for example). This pattern can be also observed in our validation site. This is an important point to consider as it means that cars will generally be grouped into one or several clusters and leave the rest of the parking empty. This property is interesting as the signal is concentrated, which makes it easier to distinguish occupied from empty areas.

Acquisition and preprocessing
We acquired Sentinel-1 GRD images from the Sentinel-Hub platform. We used Gamma0 as the backscatter coefficient, stored as floating point. Images are orthorectified, radiometrically calibrated and the thermal noise is reduced. We only keep the VH channel as it is less likely to show sinc patterns on strong reflectors.

General approach
We suppose that we have a mask P indicating where the parking is located in the image. Our images are sampled on a 10×10 meters pixels grid, so every pixel can contain up to 8 cars for typical 5×2.5 meters parking slots. We first estimate for each pixel (x, y) of an image I an occupancy ratio r(I, x, y) comprised between 0 and 1. The whole parking occupancy ratio R(I) for the image I can then be computed as: We have shown in section 2 the relationship between a pixel occupancy ratio and its brightness. We will process the problem in two ways: a simple classification where r(I, x, y) can only be 0 or 1, or a regression problem where r(I, x, y) is comprised between 0 and 1.

Occupancy estimation as a pixel classification problem
We can classify each pixel as being unoccupied -r(I, x, y) = 0 -or being occupied -r(I, x, y) = 1.
As we suppose that a pixel occupancy is positively correlated with its brightness value in the VH channel, a first approach can be to apply a simple thresholding: where v(I, x, y) is the VH value of the pixel (x, y) in the image I and t1 is a threshold value to be determined. We will name this approach C/SimpleT.
A problem with this approach is that we might count some fixed structures such as buildings or trees as occupied, since they can increase the brightness level as observed in Section 2.
Let's suppose that we are able to collect a list of acquisitions E where the parking is known to be empty or nearly empty. We are then able to compute the median M and standard deviation σ images of the empty parkings: One can take into account fixed structures by substracting M to each image I and applying a threshold on the residue: where t2 is a second threshold value to be determined. We will name this approach C/DiffT.
We can also take into account local variations in noise levels: where t3 is again a threshold value to be determined. We will name this approach C/DiffSigmaT. . Extracts of the image series in the validation site, from low occupancy to near full occupancy. The VV channel has been normalized between 0.06 and 0.18 and the VH channel has been normalized between 0.02 and 0.06. The red polygon represents the parking mask: in our case, it is the intersection of the parking shape with the Sentinel-1 image footprint as the image does not cover the whole parking. Parking occupancy was estimated by observing the pixels inside the red polygon only.
And finally, we can take into account only the general noise level: where t4 is a threshold level. This metric won't be very different from C/DiffT if only one parking is studied, but can allow to standardize results if more than one parking is monitored. We will name this approach C/DiffMedSigmaT.

Occupancy estimation as a pixel regression problem
If we assume a linear relationship between the brightness of a pixel and its occupancy ratio, we can adapt occupancy ratios shown in Section 3.3 by adding a second threshold.
where f (I, x, y) can be: As these methods are simple adaptations of classification methods, we name them: R/SimpleT, R/DiffT, R/DiffSigmaT and R/DiffMedSigmaT.

Finding threshold values
Our methods rely on threshold values that need to be specified. We propose to find those thresholds by monitoring a parking where the ground truth occupancy ratios are known and minimizing the difference between the estimated occupancy ratios and the ground truth. We develop this approach in Section 4.1.

EXPERIMENTAL VALIDATION
We shall validate and compare our proposed methods in two ways. First, we shall compare occupancy rate estimations to a validation site where the ground truth is approximately known. Secondly, we shall estimate the occupancy ratio of the 250 largest retail parkings in France between January 2018 and April 2020 and compare them to consumer and economic trends.

Validation 1
Using the camera footage taken on a validation site located in Ocean City, Maryland, we were able to estimate a fairly accurate occupancy rate. See Figure 3. There are however two limitations. First, the camera records images of the parking only every 2 minutes, so there can be up to one minute offset between a camera image and the closest Sentinel-1 acquisition. This is not a very important limitation as it is very unlikely that many cars move in one minute. A more important limitation is that the camera footage resolution is too low for annotators to discern cars that are parked far from the camera. The occupancy rate estimation therefore relies both on the number of cars counted and a visual occupancy estimation of the back of the parking. This estimation is therefore approximate.
As all methods presented in Section 3 need thresholds to be defined, we calibrated them by minimizing the loss l: where R(I) is the estimated occupancy ratio of image I and Rgt(I) is the associated ground truth occupancy ratio. Absolute difference was used instead of squared difference because first the ground truth is approximate and second we wanted to maximize the robustness, as the sample size is small. Optimization was first done through a grid search and then refined using the Nelder-Mead simplex algorithm (Nelder, Mead, 1965). We used for this effect the SciPy library (Jones et al., 2001-).
For obtaining the list of nearly empty parking images E, we selected images with a ground truth occupancy ratio less than 2%, which accounts for 10 images. This list could be easily obtained by taking into account seasonal patterns: this particular parking is near the beach, so it is generally full during summer, especially on weekends, and empty during winter.
The best thresholds and performance for our classification methods are shown in Table 1. The best obtained thresholds and performance for our regression methods are shown in Table 2. First, the C/DiffT and C/DiffMedSigmaT methods achieved the best results among the classification methods, and the R/DiffT and R/DiffMedSigmaT methods achieved the best results among the regression methods. As we were only monitoring one parking, there is no difference between both C/DiffT and C/DiffMedSigmaT, and R/DiffT and R/DiffMedSigmaT, except for the threshold values, explaining their very similar performance. The C/SimpleT and R/SimpleT methods also achieved results similar to the best methods, probably due to the fact that there is not a lot of reflective material inside the parking of the validation site. Finally, C/DiffSigmaT and R/DiffSigmaT seem to under-perform. Overall, there is not much performance difference between regression and classification methods.

Method
The relationship between the R/DiffMedSigmaT estimations and the ground truth is shown in Figure 4. To estimate the robustness of our models, we analyze how the performance varies with the threshold chosen for C/DiffMedSigmaT in Figure 5. As expected, the Mean Absolute Error is minimal on the optimal threshold estimated in Table 1. The error variance is low near the optimal threshold, which is an indication of the robustness of the model. The error remains relatively low if a higher threshold is chosen, but quickly grows when choosing lower thresholds, probably due to the SAR images noise. The occupancy rate tends to be underestimated when choosing a threshold higher than the optimal, and overestimated when choosing a lower threshold. This explains the major part of the error observed in Figure 5. However, constantly underestimating or overestimating the occupancy rate might not be problematic if the monitoring purpose is to detect trends and significant changes.
For adjusting the occupancy rate to this constant underestimation or overestimation, we can compute for each threshold t and image I the adjusted occupancy ratio R adj (I, t) from the original occupancy ratio R(I, t): R adj (I, t) = αtR(I, t), where αt is the adjustment ratio for t. αt is obtained by fitting a linear regression between R(I, t) and the associated occupancy ground truth Rgt(I) for all I (intercept set to 0). The relationship between the threshold and the adjusted Mean Absolute Error (MAE) and αt is shown in Figure 6. We observe a similar behaviour as in Figure 5, except that the MAE are overall lower, especially when choosing a lower threshold. We can also observe that the relationship between the threshold t and αt is increasing and almost linear. Overall, if we observe the range of thresholds comprised between -30% and +30% of the optimum threshold, the maximum Mean Absolute Error is 4.73% for non-adjusted estimations and 3.66% for adjusted estimations. If the general observations for the validation site stay the same for other parkings, the estimations are likely to be reliable even if the optimal threshold varies significantly between parkings.

Validation 2
The previous validation in Section 4.1 gave us a general idea of the achievable performance as well as the overall robustness of the method. However, it has two limitations. First, the number of data points is low, especially those with a high occupancy ratio (only two): although the correlation is high, it might be due to chance or other factors (weather or seasons for instance) unrelated to the presence of cars. Secondly, results were only validated on a single site and observations might not generalize well due to multiple factors (vegetation or regional specificities for example). The objective of this second validation is to check that our model is not limited to the validation site and generalizes well to other parkings.
For this purpose, we monitored 250 of the largest retail parkings in France. We didn't have access to ground truth occupancy estimations for those parkings, but checked that the observed trends are coherent with consumer and economic trends.
Parkings locations and shapes were retrieved from OpenStreet-Map (OpenStreetMap contributors, 2020). Elements with "amenity" key set to "parking" were intersected with elements with "landuse" set to "retail". The 250 parkings with the largest area were kept. The positions of all parkings are displayed in Figure 7. Although some clusters appear near Paris and other big cities, the repartition is homogeneous on most of the territory. All ascending GRD images between 2018-01-01 and 2019-12-31 were retrieved from Sentinel-hub for each parking. In total, 42539 images were collected.
In absence of known empty parking images, we started by monitoring the parking occupancy ratio of all parkings using the C/SimpleT method proposed in Section 3. We could then aggregate the relative occupancy of all the parkings using the following formula: where Ragg(d) is the aggregated occupancy ratio on day d, R(s, d) the estimated occupancy ratio of the site s on day d, and S d the set of parking sites covered by Sentinel-1 on day d. Figure 8 shows Ragg between 2018-01-01 and 2019-12-31. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition) The apparent instability of the plot is due to an important weekly pattern, as shown in Figure 9. The plot is coherent with consumer trends, as consumers tend to spend the most on Saturday and the least on Sunday (Nielsen, 2015). We estimated and verified that Sundays are the least active days. Moreover, except in December and August, most shops in France are closed on Sundays. We could therefore use this low activity to establish E, the list of empty parking images. E contains, for each parking, all Sentinel-1 images taken on nonopened Sundays. When the area had acquisitions with more than one incidence angle, we composed a set Ei for each angle and compared each image to the Ei corresponding to the image incidence angle. We could then recompute Ragg using the C/DiffT metric. See figures 10 and 11. General trends stay the same, but now the weekly pattern is much more in line with Nielsen's estimations. Figure 12 shows Ragg adjusted for weekly patterns, i.e. Ragg where each day value is divided by the corresponding weekly average shown in figure 11. A significant drop can be observed between the end of March and the end of October. This is due to the time change: Sentinel-1 images are acquired at approximately the same GMT time, but as France applies Daylight Saving Time, it means that they are acquired one hour later during France summer time, specifically between 7PM and 8PM instead of between 6PM and 7PM. Parking occupancy tends to be higher during the latter time period.
In order to adjust for this change, we compared parking occupancies 15 days before a time change compared to 15 days after the time change. Parking occupancy tends to drop 19% at the end of March whereas it tends to surge 56% at the end of October. As there are holidays at the end of October and festivities at the beginning of November, we didn't take this surge into account for our adjustment. We applied a +23.3% increase on parking occupancy estimations during summer time. Ragg adjusted for weekly patterns and summer time is shown in Figure  13.
One can observe an important increase in occupancy ratios during winter due to the festivities. We also observe significant drops on public holidays, such as the 1st of January, the 24th and 25th of December, and surges often before those holidays (for example the 23rd of December).
Most of extreme values are due to open Sundays. Since occupancy ratios are usually 3 times lower than usual during Sundays, an active Sunday creates extreme values on the weekly adjusted Ragg. Figure 14 shows the same plot but with Sundays removed. Variance is significantly reduced.
Seasonal patterns also depend on the location of monitored parkings. For instance, as shown in Figure 15, parkings on the Mediterranean coastal region -which is highly touristictend to have higher occupancy ratios during summer than during winter.
We checked the robustness of our observations by estimating Ragg for different parameters and using different methods. We show in Table 3 the absolute differences observed in Ragg adjusted for seasonality for different methods compared to the method used during this validation.
Concerning the C/DiffT method, we observe that varying the threshold by 30% only results in a maximum mean absolute error of about 0.06, which shows its robustness. Using the C/SimpleT method results in a higher difference, probably because the method doesn't use the set of empty images E.
As for C/DiffMedSigmaT, using the threshold estimated on the first validation site gives a very high error because the occupancy ratios are severely underestimated (almost nothing is detected). Difference is greatly reduced when choosing a threshold around 3.0. With an average of 0.0153 and a standard deviation of 0.0395, the median of the std images σ are overall much higher and have a very large variance compared to the one of the validation site which is around 0.0026. Most of this difference is explained because there often remains some nonempty parkings in the set of empty parkings E, as some stores are occasionally or always open on Sunday. This affects much more the std image σ than the median image M . As a consequence, the performance of this method is more dependent on the monitored parking and the quality of E than C/DiffT.
In light of the recent COVID-19 outbreak, we have plotted the occupancy ratio measured since 2020-01-01 and we observe a sharp drop from the day containment measures are put in motion. See Figure 16.

CONCLUSION
We presented in this paper a method to estimate the occupancy ratio of open car parks from Sentinel-1 image series. We have shown that there is an important correlation between estimated and ground truth occupancy ratio in a validation site. The second validation performed on 250 of the largest retail parkings in France confirms that estimations follow consumer and economic trends.
In order to get more frequent parking occupancy estimations, an axis of improvement could be to also use Sentinel-2 images. See Figure 17. In addition to issues inherent to optical satellite images, such as clouds and radiometric calibration, combining estimations from different sources -therefore with different biases and precisions -would constitute an interesting challenge.
Code and dataset can be downloaded from the following url: https://github.com/sdrdis/s1_parking_occupancy . Sentinel-1 and Sentinel-2 images of the same parking (car producer in Japan). Sentinel-1 and Sentinel-2 are acquired on the same day but at different hours, hence differences in parking occupancy are to be expected.