Juliette Blanchet
email: juliette.blanchet@univ-grenoble-alpes.fr

Jean-François Coeurjolly
email: jean-francois.coeurjolly@univ-grenoble-alpes.fr

Alexis Pellerin
email: alexis.pellerin@univ-grenoble-alpes.fr

Exploring first

Keywords: Space-time point processes, kernel estimation, Ripley's K function, global envelope test, subsampling

Introduction

Lightning strike impacts, although continuously studied by scientists in the fields of physics, climatology or statistics remains a phenomenon with an important part of randomness. Recently, there has been notable advancements in the physical investigation of the phenomenon, specifically in understanding the formation of lightning strikes. The use of the LOFAR telescope (Low Frequency Array) has played a crucial role: it enables the precise recordings of locations and times [START_REF] Hare | Lofar lightning imaging: Mapping lightning with nanosecond precision[END_REF]. We focus here on the cloud-to-ground category of lightning strikes (by opposition of the intra-cloud category when the phenomenon is concentrated only inside the lightning cloud). The goal of this paper is to analyze, from a statistical point of view, this phenomenon, thus the locations and times of impacts.

Lightning strike impacts are naturally modelled by spatio-temporal point processes, which are stochastic processes modelling events in interaction. Whether they are spatial [START_REF] Diggle | Statistical analysis of spatial and spatio-temporal point patterns[END_REF][START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF][START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF], temporal [START_REF] Daley | An introduction to the theory of point processes: volume I: elementary theory and methods[END_REF][START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF] or spatio-temporal [START_REF] Diggle | Spatio-temporal point processes: methods and applications[END_REF], point processes have known major developments these last 30 years and are now used in a large variety of fields of applications, for instance to model a disease in epidemiology [START_REF] Gabriel | stpp: an r package for plotting, simulating and analyzing spatio-temporal point patterns[END_REF], the propagation of forest fires [START_REF] Serra | Spatiotemporal log-gaussian cox processes for modelling wildfire occurrence: the case of catalonia, 1994-2008[END_REF][START_REF] Opitz | Point-process based bayesian modeling of spacetime structures of forest fire occurrences in mediterranean france[END_REF][START_REF] Raeisi | A spatio-temporal hybrid strauss hardcore point process for forest fire occurrences[END_REF], the distribution of crimes in cities (e.g. [START_REF] Mateu | Bayesian approach for modelling spatial-temporal crime data[END_REF] The dataset of interest, provided by Météorage, gives the coordinates (longitude/latitude) and the times of lightning strike impacts (in seconds) from 2011 to 2021. We focus on events occurring mainly over the French Alps which includes a part of the Italian Alps and Piemont and the Mediterranean coast, see Figure 1. The spatial observation domain corresponds to [4.43°E, 7.81°E] × [43. 1°N, 46.36°N]. As a reminder, at the equator one degree is approximately equal to 110 kilometers. We adopt a North/South division of this domain following [START_REF] Auer | Histalp-historical instrumental climatological surface time series of the greater alpine region[END_REF]. This division attempts to take into account geographical and climatological criteria. Such a division is standard in Alpine meteorological studies. It delimits the French Alps in two parts: the northern French Alps where perturbations are mainly generated by westerly flows coming from the Atlantic Ocean -we say they are under the Altlantic influence -and the southern French Alps where perturbation are mainly generated by southerly flows coming from the Mediterranean Sea -they are under the Mediterranean influence. Several studies linked extreme precipitation to the generating atmospheric influences and thus to Auer's climatological borders [START_REF] Blanchet | Retreating winter and strengthening autumn mediterranean influence on extreme precipitation in the southwestern alps over the last 60 years[END_REF] but, to the best of our knowledge, it has never been done for lightning strikes. Thus it is pertinent to see whether differences in the distribution of lightning strikes occur between Alpine regions subject to different atmospheric influences. This spatio-temporal dataset contains approximately 1.4 million of events. This dataset is illustrated in Figure 1. The red curve represents the Auer's climatological border between regions called 'North Alps' and 'South Alps'.

Locations of impacts are superimposed on the altitude map for a better visualization.

Times, locations as well as the number of total events of lightning strikes are random and are therefore modelled by a spatio-temporal point process. In this paper, we intend to understand on the one hand the inhomogeneity of events across time and space and on the other hand dependence between these events. Usually, the first question is addressed by analyzing the (first-order) temporal, spatial or spatio-temporal intensity function and it is pertinent to investigate whether this intensity is separable in space and time (see e.g. [START_REF] Møller | Aspects of second-order analysis of structured inhomogeneous spatio-temporal point processes[END_REF]. The second question can be tackled using higher-order summary statistics like the Ripley's K function or the J function (see e.g. [START_REF] Diggle | Statistical analysis of spatial and spatio-temporal point patterns[END_REF][START_REF] Cronie | A j-function for inhomogeneous spatio-temporal point processes[END_REF] that measure departures from the homogeneous or inhomogeneous Poisson point process, the reference process modeling independent events in time and space. All these methodologies require fine tuning of hyperparameters (like bandwidth for kernel type estimators) or quite a large number of simulations for tests based on Monte-Carlo replications and global envelopes (e.g. [START_REF] Myllymäki | Global envelope tests for spatial processes[END_REF][START_REF] Ghorbani | Testing the first-order separability hypothesis for spatio-temporal point patterns[END_REF]. Due to the considerable volume of events and with the aim of minimizing computational time, this paper focuses on exploring these two questions through the use of subsampling, a technique that has recently regained significance in spatial point pattern analysis [START_REF] Cronie | A cross-validation-based statistical theory for point processes[END_REF].

In the rest of this paper, we first present in Section 2 a general background and notation on spatio-temporal point processes. Moments and summary characteristics for spatio-temporal point processes and versions of these characteristics after subsampling are presented. Section 3 presents results on nonparametric intensity estimation, first-order separability tests and estimation of space-time subsampled Ripley's K functions. Finally, a conclusion is presented in Section 4

while Appendix A illustrates an idea presented in Section 2.4, which consists in subsampling an inhomogeneous point process to make it homogeneous.

2. Spatio-temporal point processes and subsampling

Background on point processes

We view a point process on a (complete metric) space set, say S, as a locally finite configuration of events in S (see e.g. [START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF]; [START_REF] Daley | An introduction to the theory of point processes: volume I: elementary theory and methods[END_REF]). When

S = R + (resp. S = R d , S = R d × R +
) the stochastic models are often referred as temporal (resp. spatial, spatio-temporal) point processes. We focus in this paper on planar spatio-temporal point processes X st defined on S = R 2 × R + . Our observation domain is denoted by W × T where W ⊂ R 2 is a bounded spatial domain (here the French Alps) and T ⊂ R + the period of observation (here 2011-2021).

Therefore, the data consists in a set x st = {y i = (x i , t i ), i = 1, . . . , n} where x i ∈ W and t i ∈ T respectively correspond to the observed location and time of the ith lightning strike.

Thus y i stands for the space-time location. Without loss of generality, we order events by their observed times, t 1 < ... < t n . Locations, times and n are realizations of randoms variables and we assume that X st is simple, meaning that (x i , t i ) ̸ = (x j , t j ) for any i ̸ = j.

We let X t and X s denote the point process aggregated in time or space, that is

X t = {t : y = (x, t) ∈ X st , x ∈ W } and X s = {x : y = (x, t) ∈ X st , t ∈ T } . We let N (A × B) = # {(x, t) ∈ X st ∩ (A × B)} denote the counting variables, that is number of events in A × B where A × B ⊂ R 2 × R + .
With a slight abuse of notation, we also denote

N (A) and N (B) by N (A) = # {x ∈ X s ∩ A} and N (B) = # {t ∈ X t ∩ B}.
A point process is stationary (resp. isotropic) if its distribution is invariant under translations (resp. rotations).

The distribution of X st (as well as those of X t , X s ) is characterized by the void probabilities, the finite dimensional distributions of counting variables or via their generating moment functionals (see e.g. [START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF]). These are usually difficult to establish. Intensity functions and summary statistics are simpler to define and estimate. Some of them are defined in the next two sections.

First-order intensity functions

Assume that the processes X st , X s and X t have first-order intensities denoted by λ st , λ s and λ t . These functions can be interpreted for any y = (x, t) with x ∈ R 2 and t ∈ R + by

λ st (y) = lim |dy|→0 E{N (dy)} |dy| ; λ s (x) = lim |dx|→0 E{N (dx)} |dx| ; λ t (t) = lim |dt|→0 E{N (dt)] |dt| (1)
where dy = dx × dt. By application of the Campbell theorem (see [START_REF] Daley | An introduction to the theory of point processes: volume I: elementary theory and methods[END_REF]) these intensities are linked by

λ s (x) = T λ st (x, t)dt and λ t (t) = W λ st (x, t)dx (2)
and are obviously linked to counting variables by

E{N (A × B)} = A×B λ st (y)dy, E{N (A)} = A λ s (x)dx, E{N (B)} = B λ t (t)dt (3)
for any bounded Borel set A × B ⊂ R 2 × R + . A spatio-temporal point process is said to be homogeneous if λ st (y) = λ for any y ∈ R 2 × R + . Otherwise, the process is said to be inhomogeneous. In this case, the intensity could depend only on x, only on t or both on x and t. And in the latter several authors (see e.g. [START_REF] Møller | Aspects of second-order analysis of structured inhomogeneous spatio-temporal point processes[END_REF]) investigate a separability hypothesis of the intensity function. A spatio-temporal point process is said to be first-order separable if λ st (x, t) can be factorized as

λ st (x, t) = λ sep st (x, t) = λ 1 (x)λ 2 (t) for any (x, t) ∈ W × T (4)
where λ 1 and λ 2 are non-negative and measurable functions respectively on R 2 and R + . Under the hypothesis of first-order separability, i.e. if (4) holds, the intensity of X s and X t can respectively be obtained as

λ s (x) = λ 1 (x) T λ 2 (t)dt and λ t (t) = λ 2 (t) W λ 1 (x)dx.
(5) Thus, by combining (4) and ( 5) under the first-order separability assumption the spatio-temporal intensity writes

λ sep st (x, t) = λ s (x)λ t (t) ν with ν = E{N (W × T } = W ×T λ st (x, t)dxdt.
Following [START_REF] Ghorbani | Testing the first-order separability hypothesis for spatio-temporal point patterns[END_REF], we introduce the natural summary functions S st (x, t) = 100 λ st (x, t)/λ sep st (x, t) = ν λ st (x, t)/{λ s (x)λ t (t)} and their spatial and temporal averages Estimates of λ t , λ s , λ st , S t , S s and S st and tests of first-order separability using global envelopes [START_REF] Myllymäki | Global envelope tests for spatial processes[END_REF] are investigated in Section 3.

S t (t) = |W | -1 W S st (x, t)dx and S x (x) = |T | -1 T S st (x,

Higher-order intensity functions and summary statistics

If the kth (k ≥ 1) moment measure of X st is absolutely continuous with respect to the Lebesgue measure of R 2 × R + , the kth order intensity exists and is defined as

λ (k) st (
st (y 1 , . . . , y k ) = λ st (y 1 ) . . . λ st (y k ). The pair correlation is an index of departure to Poisson assumption and is defined for any distinct [START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF])

y 1 , y 2 ∈ R 2 × R + by (see e.g.
g st (y 1 , y 2 ) = λ (2) st (y 1 , y 2 ) λ st (y 1 )λ st (y 2 )
where we use the convention a/0 = 0 for any a ≥ 0. If X st is stationary (resp. isotropic) then λ st is constant over space and time and g st depends only on y 2 -y 1 (resp. ∥y 2 -y 1 ∥). For inhomogeneous (thus non stationary) point processes and to be closer to the space-time nature of X st several frameworks and assumptions exist. Among them, we consider [H 2 ] and [H

′ k ] [H 2 ] λ (2)
st exists, g st is invariant under translations and there exists a function ḡ such that

g st (y 1 , y 2 ) = g(y 2 -y 1 ) = ḡ(∥x 2 -x 1 ∥, |t 2 -t 1 |). [H ′ k ] λ (k)
st (y 1 , . . . , y k ) exists and satisfies for any pairwise distinct y 1 , . . . ,

y k ∈ R 2 × R + and any a ∈ R 2 × R + λ (k) st (y 1 , . . . , y k ) λ st (y 1 ) . . . λ st (y k ) = λ (k) st (y 1 + a, . . . , y k + a) λ st (y 1 + a) . . . λ st (y k + a)
.

[H 2 ] assumes that the pair correlation depends only on the spatial and temporal distances between events. It is particularly well-suited for spatio-temporal point processes where time and space typically serve distinct roles and operate on different scales. Under [H 2 ], Gabriel et al.

(2013) introduced the space-time Ripley's K function as a very natural extension of the stationary version by

K inh,st (r, τ ) = Br,τ ḡ(∥x∥, |t|)dxdt (7)
where B r,τ is the cylindrical ball centered at 0 given by B r,

τ = {y = (x, t) ∈ R 2 × R + : ∥x∥ ≤ r, |t| ≤ τ }. Under the Poisson case, g st = ḡst = 1 and K inh,st (r, τ ) = |B r,τ | = 2τ πr 2 .
Inhomogeneous versions of standard statistics such as the empty space function F , the nearest-neighbour distribution function G and the J function are less straightforward to derive. [START_REF] Cronie | A j-function for inhomogeneous spatio-temporal point processes[END_REF] consider intensity reweighted moment stationary (IRMS) spatio-temporal point process, that is models satisfying [H ′ k ] for any k ≥ 1, and such that λ > 0 and lim sup

k→∞ λk k! I k (r, τ ) 1/k < 1 (8)
where λ = min y∈R 2 ×R + λ st (y) > 0 and where for any k ≥ 1

I k (r, τ ) = Br,τ • • • Br,τ λ (k) st (y 1 , . . . , y k ) λ st (y 1 ) . . . λ st (y k ) dy 1 . . . dy k Ĩk (r, τ ) = Br,τ • • • Br,τ λ (k+1) st (0, y 2 , . . . , y k+1 ) λ st (0) . . . λ st (y k+1 ) dy 2 . . . dy k+1 .
Under this framework, [START_REF] Cronie | A j-function for inhomogeneous spatio-temporal point processes[END_REF] use series expansions of the moment generating functionals to propose the following extensions of F, G, J, provided here with little details, for any r, τ > 0

1 -F inh,st (r, τ ) = 1 + k≥1 (-λ) k k! I k (r, τ ) 1 -G inh,st (r, τ ) = 1 + k≥1 (-λ) k k! Ĩk (r, τ ) J inh,st = 1 -G inh,st (r, τ ) 1 -F inh,st (r, τ ) .
To see the intuition, it can be shown (see again [START_REF] Cronie | A j-function for inhomogeneous spatio-temporal point processes[END_REF]) that for stationary spatio-temporal point process, 1 -F inh,st (r, τ ) = P(N (B r,τ ) > 0) and 1 -G inh,st = P(N (B r,τ ) > 0 | 0 ∈ X st ) which correspond to standard definitions of these functions. Finally, under the Poisson assumption, we may check that 1 -F inh,st (r, τ ) = 1 -G inh,st (r, τ ) = exp(-λ|B r,τ |) and J inh,st (r, τ ) = 1. In the rest of the paper, we assume that X st satisfies [H 2 ],

[H ′ k ] for any k ≥ 1 and (8).

Subsampling point processes

The nonparametric estimation of intensity functions, tests of first-order separability, goodnessof-fit tests of inhomogeneous Poisson point processes require a tuning of few parameters (like bandwidths parameters) and/or simulations which can lead to severe computational cost. This mainly comes from the highly inhomogeneous spatio-temporal nature of the lightning strikes dataset (there are many areas and large periods of time where no lightning strike is observed)

and the large number of observed events. To reduce these costs, it does therefore make sense to consider subsamples of the dataset.

Subsampling point processes has regained popularity recently (see e. we define the subsampled version X π st as the thinning with probability 1 -π (or conversely with retaining probability π) as

X π st = { y ∈ X st : ε(y) = 1 } .
When π(•) = π 0 is constant over space and time, X π0 st is a space-time independent subsampling of X st . In this section we briefly review properties of X π st . When we add the superscript π to a characteristic, for instance λ (k),π st , g π st we mean the kth order intensity function, the pair correlation function of X π st . Iterated versions of Campbell theorem allow to prove, see e.g. [START_REF] Cronie | A cross-validation-based statistical theory for point processes[END_REF], that

λ (k),π (y 1 , . . . , y k ) = π(y 1 ) . . . π(y k )λ (k) (y 1 , . . . , y k ) = π k 0 λ (k) (y 1 , . . . , y k ) (9)
where the latter holds if π(•) = π 0 . From (9), we deduce the following proposition.

Proposition 1. Under the previous assumptions on X st and general notation, we have the following statements (i) X π st necessarily satisfies the same assumptions as X st , namely

[H 2 ], [H ′ k ] for any k and (8) if π = inf y π(y) > 0. (ii) If π(y) = π 0 , S π0 st (x, t) = S st (x, t), S π0 s (x) = S s (x) and S π0 t (t) = S t (t).
(iii) The pair correlation and the space-time Ripley's K function satisfy for any y 1 ̸ = y 2 and any r, τ > 0

g π st (y 1 , y 2 ) = g st (y 1 , y 2 ) and K π inh,st (r, τ ) = K inh,st (r, τ ). (iv) From (i), X π st is an IRMS, for any r, τ > 0, I π k (r, τ ) = I k (r, τ ), Ĩπ k (r, τ ) = Ĩk (r, τ ) which yields 1 -F π inh,st (r, τ ) = 1 + k≥1 (-λπ) k k! I k (r, τ ) (10) 1 -G π inh,st (r, τ ) = 1 + k≥1 (-λπ) k k! Ĩk (r, τ ). (11)
Proof of Proposition 1 is quite straightforward and essentially follows from (9). Note that (ii) is the only result for which the subsampling is homogeneous. This result would not be true even if the retaining probability field is separable in space and time, that is if π(y) = π(x, t) = π s (x)π t (t).

Results ( 9), (i)-(iii) show that we can recover characteristics like intensity functions of X st from the ones of X π st . So to estimate g st or K inh,st , we simply estimate g π st and K π inh,st . To estimate λ st (y), we can estimate λ π st and set λst (y) = λπ st (y)/π(y). Result (iv) was proved by [START_REF] Cronie | A j-function for inhomogeneous spatio-temporal point processes[END_REF]. It tells that the subsampling has a more complex effect on distance-based summary statistics like F inh,st , G inh,st and J inh,st . This remark would still apply even if X st were stationary and/or if we consider a homogeneous subsampling.

Going back to the general objective of this section which is to reduce computational cost of nonparametric estimation, we propose two subsampling strategies:

(a) π(y) = π 0 with π 0 ∈ (0, 1). In particular, we view π 0 as a small positive real number.

(b) Let µ > 0 and π(y) = µ λst(y) 1(y ∈ ∆ µ ) where ∆ µ = {y ∈ R 2 × R + : λ st (y) ≥ µ}.
Note that if one sets the expected number of points ν π of X π st , we can set π 0 and µ as follows

ν π = E{N π (W × T )} = W ×T π(y)λ st (y)dy =    π 0 E{N (W × T )} (a) µ|∆ µ | (b)
For (a), E{N (W × T )} is estimated by n the observed number of data points and so π 0 = ν π /n.

For (b), the problem may not have a solution (or a unique solution) but µ could be obtained using an optimization procedure (see Appendix A for more details).

Strategy (a) is very natural. We use no information from X st and subsample independently of space and time. Strategy (b) is more tricky. Indeed, in particular

λ π st (y) = π(y)λ st (y) = µ1(y ∈ ∆ µ ).
So X π st | ∆µ becomes a homogeneous point process with intensity µ on ∆ µ . This strategy suffers from the obvious drawback that λ st is unknown and we believe this strategy is uninteresting to estimate λ st or to test the first-order separability. However, for the estimation of K inh,st (r, τ ), the interest is clear: X π st is homogeneous on ∆ µ and has the same second-order structure as X st . So, for instance, to construct global envelopes of K π inh,st under the Poisson assumption,

200
we simply have to simulate homogeneous Poisson point processes on ∆ µ . Appendix A explores these ideas: we propose a practical procedure and simulation study (only for planar spatial point processes). We show that strategy (b) has some strong merit and deserves to be investigated in a future research. However, we did not pursue this idea for the lightning strikes dataset. The main difficulties are that the intensity λ st is highly inhomogeneous (lots of peaks and empty 'space-time' areas). This makes the set ∆ µ really non convex and thus the estimation of the space-time (homogeneous) Ripley's K function really complex and unstable as it must take into account border correction.

We end this section with the following proposition which shows the impact of strategy (a), with a small value of π 0 , for the estimation of distance-based summary statistics. This result provides a more formal approximation of Cronie and Van Lieshout (2015, bottom of p.16).

Proposition 2. Let X π0 st be a subsampled version of X st with π(•) = π 0 > 0. We assume [H 2 ], [H ′ k ] for any k ≥ 1 and (8). Then under the notation of Proposition 1, for any r, τ > 0, we have the following statements as π 0 → 0

(i) 1 -F π0 inh,st (r, τ ) = 1 -F π0,Poisson inh,st (r, τ ) 1 + π 2 0 λ2 2 I 2 (r, τ ) + o(π 2 0 ) (12) (ii) 1 -G π0 inh,st (r, τ ) = 1 -G π0,Poisson inh,st (r, τ ) 1 + λ2 π 2 0 2 Ĩ2 (r, τ ) + o(π 2 0 ) (13) (iii) J π0 inh,st (r, τ ) = 1 -λπ 0 {K inh,st (r, τ ) -|B r,τ |} + λ2 2 π 2 0 (K inh,st (r, τ ) -|B r,τ |) 2 + Ĩ2 (r, τ ) -I 2 (r, τ ) + o(π 2 0 ) (14)
Proof. Let (u n ) n be a sequence of real numbers such that u n = o(1/n) as n → ∞, we leave the reader to prove that for any c > 0

1 + c n + u n n = exp(c) 1 + 1 2 nu 2 n + o(1/n) . (15)
Now, consider (10). Since from ( 8)

k≥3 (-λπ 0 ) k k! I k (r, τ ) ≤ π 3 0 k≥ (-λπ 0 ) k k! I k (r, τ ) = O(π 3 0 ) we deduce that 1 -F π0 inh,st (r, τ ) = 1 -λπ 0 I 1 (r, τ ) + λ2 2 π 2 0 I 2 (r, τ ) + o(π 2 0 ).
Using (15) we then obtain, by noticing that

I 1 (r, τ ) = |B r,τ | 1 -F π0 inh,st (r, τ ) 1/π0 = exp -λ|B r,τ | 1 + λ2 2 π 0 I 2 (r, τ ) + o(π 0 )
which leads to (i) by reminding of course that 1 -F π0,Poisson inh,st (r, τ ) = exp(-π 0 λ|B r,τ |). We now consider (11). Since (8) also implies that k≥3

(-λπ 0 ) k k! Ĩk (r, τ ) = O(π 3 0 )
we have using ( 15)

1 -G π0 inh,st (r, τ ) 1/π0 = exp -λ Ĩ1 (r, τ ) 1 + λπ 2 0 2 Ĩ2 (r, τ ) + o(π 2 0 ) 1 -G π0 inh,st (r, τ ) 1/π0 = 1 -G π0,Poisson inh,st exp -λπ 0 ( Ĩ1 (r, τ ) -I 1 (r, τ )) × 1 + λ2 π 2 0 2 Ĩ2 (r, τ ) + o(π 2 0 ) .
The result is obtained by noticing that Ĩ1 (r, τ ) = Br,τ g st (y)dy = K inh,st (r, τ ) and by applying a Taylor expansion of the exponential function. Finally, (14) follows easily from ( 12)-( 13).

Thus, when π 0 is small, 1-F π0 inh,st (r, τ ) = 1-F π0 inh,st (r, τ )(1+o(π 0 )) = (1-exp(-λπ 0 B r,τ ))(1+ o(π 0 )) and J π0 inh,st = 1 + O(π 0 ), so X π0 st is hard to distinguish from an inhomogeneous Poisson point process. One could be tempted to define Jπ0 inh,st (r, τ ) = (1 -J π0 inh,st (r, τ ))/π 0 to discard the first-order term but in that case

Jπ0 inh,st (r, τ ) = λ K π0 inh,st (r, τ ) -|B r,τ | + o(π 0 ).
Hence, upto o(π 0 ), Jπ0 inh,st (r, τ ) only measures the difference between the inhomogeneous K function and the one under the inhomogeneous Poisson assumption. Proposition 2 combined with the fact that λ is extremely small on the lightning strikes dataset making the estimation of J inh,st very unstable, we have decided to skip distance-based summary statistics in the data analysis.

To summarize this section, among different strategies for subsampling, the one assuming no information a priori seems the more appropriate to estimate characteristics of X st . We typically use π 0 ≈ 2.5%, so that the number of points in X π0 st is close to 35,000 points, which is still quite large. According to Propositions 1-2, Section 3 presents empirical results for the estimation of λ st , S st , S s , S t and K inh,st .

Results on the dataset

Intensity estimation

In this section, we present empirical results for the nonparametric estimations of λ t , λ s and λ st . We consider the kernel intensity estimator with edge effect correction (see e.g. [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF]). To estimate λ st instead of a 3-dimensional kernel, we follow [START_REF] Ghorbani | Testing the first-order separability hypothesis for spatio-temporal point patterns[END_REF] and use a product of two kernels (one in space and one in time). For any y = (x, t) ∈ W × T , we define for two bandwidths parameters b

s , b t λst (y; b s , b t ) = y ′ =(x ′ ,t ′ )∈Xst k s,bs (x ′ -x) e s (x ′ ; b s ) k t,bt (t ′ -t) e t (t ′ ; b t )
where

k s,bs = b -2 s k s (•/b s ), k t,bt = b -1 t k t (•/b t ) and
where k s (resp. k t ) is a two-dimensional (resp. one-dimensional) kernel function. The terms e s and e t serve as edge correction factor and we use the Diggle's correction [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF] given by

e s (x ′ ; b s ) = W k s,bs (x -x ′ )dx and e t (t ′ ; b t ) = T k t,bt (t -t ′ )dt.
Estimates of λ s and λ t can be obtained from λst or directly from X s and X t . We use the latter and define for any y = (x, t) ∈ W × T and two bandwidths parameters bs and bt λs (x; bs ) =

x ′ ∈Xs k s, bs (x ′ -x) e s (x ′ ; bs ) and λt (t; bt ) =

t ′ ∈Xt k t, bt (t ′ -t)
e t (t ′ ; bt ) .

As illustrated by [START_REF] Bivand | Applied spatial data analysis with R[END_REF], the choice of the kernel to be used has less impact than the smoothing parameters. We simply used isotropic Gaussian kernels: thus k s and k t depend only on ||x ′ -x|| and |t ′ -t| respectively. We select the smoothing parameter bt as proposed by [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF]. Regarding the choice of bs , we follow the simulation done by [START_REF] Cronie | A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions[END_REF] and use a similar approach. They suggest to select bs as the parameter achieving the smallest squared inverse residuals value (see [START_REF] Baddeley | Residual analysis for spatial point processes (with discussion)[END_REF]) given by

L( bs ) = x ′ ∈Xs 1 λs (x ′ ; bs ) -|W | 2 . ( 16
)
In this procedure, λs (•; bs ) has to be estimated precisely at data points for any bs . This can be extremely computationally expensive if one uses the estimator proposed by [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF] for that task. Instead, we use a subsampling strategy combined with a 10-fold cross-validation.

We first substitute X s with a subsampled version X π s (obtained with π(•) = π 0 = 2.5%). Then, in ( 16) we estimate λ s using 90% of data from X π s and evaluate the loss function at the 10% remaining points. We repeat the estimation of the loss function 10 times, then compute the average loss and finally the optimal bs . We also repeat the subsampling procedure 50 times and average the bs .

Estimates of λ t and λ s are evaluated over a regular grid of respectively 1,000 times points between 01/01/2011 -00:00:00 and 31/12/2021 -23:59:59 and a 256 × 256 pixels grid. The estimation of λ st is more expensive and we use a homogeneous subsampling with π 0 = 2.5% (such that ν π ≈ 35, 000) and Proposition 1. Also, to reduce computational complexity, we use b s = bs and b t = bt for the estimation of λ st . Figures 23depict empirical results for estimates of λ s and λ t only.

Let us comment Figure 2 which explores variations between seasons and between the northern and southern parts of the French Alps. The impact rate in the northern Alps seems to be, on average, more than twice as low as that in the southern Alps, across all years and seasons. This is probably related to the Mediterranean influence that is predominant in the southern region [START_REF] Blanchet | Retreating winter and strengthening autumn mediterranean influence on extreme precipitation in the southwestern alps over the last 60 years[END_REF], a region known to generate more convective events, which are associated to electrical activity. Furthermore, there is a significant year-to-year variability in the estimated intensity. Specifically, 2017 appears to be an exceptional year characterized by a low number of lightning strikes throughout the year, with an average intensity about four times lower than that observed in the three areas. The year 2017 was indeed characterised by predominant anticyclonic conditions responsible for a significant shortfall in rainfall over almost the entire country. Conversely, 2021 recorded nearly three times the average number of impacts in the three areas. In 2021, lightning strikes were more frequent in the northern regions during spring and in the southern regions during autumn. Considering that these regions experience distinct atmospheric influences (Atlantic vs. Mediterranean), this highlights the role of atmospheric circulation in the spatiotemporal variability of lightning strikes. Overall, the impact intensity is lower in winter and higher in summer (with summer intensity nearly 40 times the average winter intensity), attributed to the predominance of convective events during the summer season. We now comment Figure 3 which displays estimates of λ s (x) over a 256 × 256 pixel grid. We estimate λ s for the whole period as well as every year from 2011 to 2021. The Italian region seems to be prone to a significantly higher number of lightning strike impacts, with an intensity exceeding ten times the regional average.The Piemont region is indeed recognized for experiencing "East returns," which refer to perturbed weather patterns leading to intense storms and the potential for significant precipitation, as highlighted in the work of (Blanchet et al., 2021a). In general, Auer's border effectively delineates the spatial intensity of lightning strikes, revealing higher intensities in the southern part of the Alps compared to the northern part. This is in agreement with the conclusions obtained from Figure 2. Moreover, the annual maps indicate significant year-to-year variability, particularly in 2017, which stands out as a markedly different year with an intensity nearly six times lower than the overall country's average, as depicted in Figure 2.

First-order separability testing

Based on the previous section, we now investigate the first-order separability hypothesis of the intensity function, that is the null hypothesis given by (4). We use the summary statistics S st , S s and S t presented in Section 2, see ( 6). Since we estimate λ st using subsampling, we proceed similarly for these summary statistics, as guaranteed by Proposition 1. Then, we follow the statistical procedure proposed by [START_REF] Ghorbani | Testing the first-order separability hypothesis for spatio-temporal point patterns[END_REF] based on permutations to test the null hypothesis and in particular to have replications of S st , S t and S s under the null. Finally, we use global envelope tests for S t and S s [START_REF] Myllymäki | Global envelope tests for spatial processes[END_REF] to have a Monte-Carlo valid p-value and graphical interpretation of departures to the null. We have reproduced the whole 300 of this procedure on two versions of X π0 st to investigate the robustness of findings. Figure 4 presents empirical results. Even if there are some (expected) slight differences between the two versions of X π0 st , the conclusions are clear. In the two versions of X π0 st , there is a strong evidence that the spatio-temporal intensity function is not separable in space and time. It can be also be pointed out that the first-order separability test is rejected mainly during summers periods mainly either in the southern Alps or in the part of the northern French Alps that experiences East return events (the Savoie and Haute-Savoie departments close to the Italian border). This may be attributed in both cases to the occurrence of convective storms showing strong electrical activity over short time periods and limited areas, inducing strong space-time dependencies. Figure 4: Results for the first-order separability testing procedure obtained from two subsampled versions of X π 0 st (first and second rows). 95% global envelopes tests using the ERL procedure were used where permutations are used to generate data under the null. The considered summary statistics are St(•) (left, solid curve) and Ss(•) (right, latent image). Envelopes are depicted only for St and red dots indicate times or pixels for which the estimated summary statistics is outside the envelopes. Light blue rectangles corresponds to summer seasons (from June to end of September). The adjusted global p-value for each of these envelope test is p = 0.05%.

Space-time Ripley's K function

We now dig in the higher-order structure of X st . As specified at several places, we consider only subsampled versions of X π0 st with π 0 = 2.5% (so that ν π ≈ 35, 000). As a consequence of Section 2.4 and in particular in light of Propositions 1-2, we consider only the space-time inhomogeneous Ripley's K as relevant summary statistic, that is K inh,st (r, τ ) = K π0 inh,st (r, τ ) given by ( 7). Let r ∈ [r m , r M ] and τ ∈ [τ m , τ M ]. To derive a simple graphical interpretation, we also suggest to consider and estimate the following one-dimensional functions K t (τ ) and K s (r)

given by

K inh,t (τ ) = 1 3π(r 3 M -r 3 m ) W K inh,st (r, τ )dr (17) K inh,s (r) = 1 τ 2 M -τ 2 m T K inh,st (r, τ )dτ. ( 18
)
The idea behind ( 17)-( 18) is that under the Poisson assumption K inh,t (τ ) = 2τ and K inh,s (r) = πr 2 which respectively correspond to the Ripley's K function of a one-dimensional and twodimensional inhomogeneous Poisson point process. Note that K inh,t (τ ) and K inh,s (r) do not correspond to the Ripley's K functions of X t and X s . Estimates of K inh,t (τ ) = K π0 inh,t (τ ) and K inh,s (r) = K π0 inh,s (r) ensue from the estimate of K inh,st (r, τ ). As seen from Proposition 1, we use a (homogeneous) subsampled version X π0 st . The estimator, proposed by [START_REF] Gabriel | stpp: an r package for plotting, simulating and analyzing spatio-temporal point patterns[END_REF], is given for some r, τ > 0 by

Kinh,st (r, τ ) = 1 |W × T | y,y ′ ∈X π 0 st e(y, y ′ ) λs,t (y) λs,t (y ′ ) 1(y -y ′ ∈ B r,τ ).
where e(•, •) is an edge-correction factor (see e.g. [START_REF] Møller | Statistical inference and simulation for spatial point processes[END_REF]; [START_REF] Gabriel | Estimating second-order characteristics of inhomogeneous spatio-temporal point processes: influence of edge correction methods and intensity estimates[END_REF]) and where λst is obviously a preliminary estimation of the spatio-temporal intensity (eventually itself estimated using a subsampled version of X st ).

Then, from the estimation of K inh,t and K inh,s we construct a combined global envelope test (using the extreme rank length procedure) proposed by [START_REF] Myllymäki | Global envelope tests for spatial processes[END_REF]. To construct envelope, we generate B = 199 simulations of inhomogeneous Poisson point processes on W × T .

We evaluate Kinh,t (τ ) (resp. Kinh,s (r)) for 50 values of τ (resp. r) from τ m = 0 (resp. r m = 0)

to τ M = 0.75% × |T | (resp. r M = 20% × |W | 1/2
). The choice of r M follows from rule of thumb for spatial point processes [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF]. The value of τ M represents a one-month period approximately. Estimates Kinh,t and Kinh,s and envelopes are represented in Figure 5.

To facilitate the interpretation x-axis for Figure 5 (left) we have rescaled time events such that all events occur in time in the interval (0, 1).

Conclusion

In this paper, we estimate first and second order structure of the inhomogeneous spatiotemporal point process given by lightning strike impacts. We considered subsampling strategy to tune kernel estimation, to estimate the null distribution of the first-order separability test, to estimate the inhomogeneous K function and its behavior under the Poisson assumption. The results are crystal clear, this point pattern is highly inhomogeneous both in time and space.

This dependence is time and space is complex and the process seems to exhibit a high degree of clustering in time and space. These conclusions were addressed using largely reduced subsampled which allows us to reduce significantly computational time and cost. By passing we have shown that a small value of the constant retaining probability π 0 prevents from understanding distancebased summary statistics which are almost indistinguishable from an inhomogeneous Poisson point process. This research leads to many perspectives. From a theoretical point of view, some are presented in Appendix A and concern the strategy (b) to subsample (ie. transform an inhomogeneous point process into a homogeneous one). Others could concern the precision of subsampled estimators like λπ st or Kπ inh,st . Propositions 1-2 have shown that subsampling methodology is a correct procedure but understanding the rate of convergence of nonparametric estimators with the retaining probability function could be of great importance.

From a practical point of view, after the exploratory analysis presented in this paper, the natural next step is to model more specifically the intensity function. We have already presented a few spatial covariates like the altitude map, the binary map induced by the Auer division between the North Alps and South Alps. We also have at our disposal several other spatiotemporal covariates like temperature, levels of pressure, heights of precipitation, etc. These extra information consitute more than 25Go of data and raise many other numerical and practical challenge we expect to investigate in a future research.

which is a valid test in this context. Poisson point processes X s are generated in W = [0, 1] 2 with intensity function from the two following models

λ s (x) = β    exp -ϕ -1 (1 + |x 1 -1/2| × |x 2 -1/2|) (cross) exp ϕ -1 × (| sin(πx 1 ) + | sin(πx 2 ) (sin)
where x = (x 1 , x 2 ) ∈ W , ϕ = 0.02 (resp. 0.2) for the cross model (resp. sin model) and where we set β such that ν = EN (W ) = 50, 000 points. Finally we conduct the simulation such that ν π = 500 so we keep on average 10% from the initial pattern. it can be checked that for these models and simulations the losses functions equal 0 at these minima. Figures A.8-A.9 present simulated patterns. First rows correspond to simulations of X s on W = [0, 1] 2 with on average 50,000 points. Second rows are subsampled version X π s obtained with the strategy (b) where we used the theoretical retaining probability π and thus the true set ∆ µ . Thus the grey areas correspond to W \ ∆ µ . Finally, third rows correspond to simulations of X π s obtained with the data-driven procedure described above where we both estimate π and ∆ µ . Visually speaking, we observe that the subsampled versions X π s and X π s tend to produce homogeneous patterns respectively on ∆ µ and ∆µ . To confirm this, we also perform the quadrat test procedure [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF] (where the number of tiles follows from the rule of thumb from the spatstat package) for all point patterns (so for X s , X π s and X π s ). The title of each graph contains the observed number of points and the p-value of this test (accurate to 0.01%).

Obviously, the rest rejects (with high confidence) the homogeneous Poisson point process (on W ) assumption for the patterns X s (p-values are all smaller than 0.01%). Values show that the null hypothesis (homogeneous Poisson on ∆ µ or ∆µ ) cannot be rejected at level 5% for simulations of X π s and X π s . This simulation study suggests that the proposed algorithm which aims to subsample an inhomogeneous point process to make it homogeneous seems possible. It does raise many questions. Among them : (a) are there theoretical guarantees that the estimated loss function has a unique minimum, that the value of the loss function at this minimum is zero? (b) are the findings on these models reproducible for a larger class of models and in particular for spatio-temporal guarantees? (c) can we quantify the gain of having a homogeneous subsample compared to an inhomogeneous one? Definitely, these interesting questions will be the topic of further research. Vertical dashed lines correspond to minima of these functions.

Figure 1 :

 1 Figure 1: Locations of lightning strike impacts aggregated per day from June 2nd to June 7th 2018 (a quite active week) observed in the French Alps study domain, which includes a part of the Mediterrannean and Italian Alps.

  g.[START_REF] Chiu | Stochastic geometry and its applications[END_REF];[START_REF] Cronie | A cross-validation-based statistical theory for point processes[END_REF]) in particular in the context of statistical learning, cross-validation technique, variance estimation, etc. Independent subsampling corresponds to the process of thinning (or conversely retaining) a point from an initial point pattern. Let π : R 2 × R + → [0, 1] and (ε(y), y ∈ R 2 × R + ) a random field of independent Bernoulli distributions with parameter π(y),

Figure 2 :Figure 3 :

 23 Figure 2: Nonparametric temporal intensity estimated by study area (Northern Alps, Southern Alps, All the French Alps) for the 2011-2021 period; Nonparametric temporal intensity estimated by study area (Northern French Alps, Southern French Alps, All the French Alps) by season (Winter=December to March, Spring=April to May, Summer=June to August, Autumn=September to November).

Figure 5 :

 5 Figure5: Combined global envelope test based on the summary functions K inh,t (τ ) and K inh,s (r). We have rescaled time events in (0, 1) to have an easier interpretation of the x-axis for K inh,t (τ ) (in particular 0.075 corresponds approximately to a one-month period). Observed curves have been truncated on the y-axis for a better visualization. Insets represent the non-truncated versions.

  Figure A.6 illustrates the intensity function models. Figures A.8-A.9 present empirical results for two simulated patterns x s from these models (point patterns are illustrated in the first rows of these figures). For each simulation, we adopt the previously described algorithm. Figure A.7 shows that the loss functions L(µ) for both models are well estimated. The minimum of the loss function is also well estimated. And

  Figure A.6: Intensity functions of the two models considered in the simulation. The cross model (left) concentrates points around the x-axis or the y-axis, while the sin model tends to produce patterns in four clusters. Colors are in log-scale.

Figure A. 8 :

 8 Figure A.8: First row present two simulations with on average 50,000 points of the cross model. Second row (resp. third row) present two subsampled versions X π s with the theoretical retaining probability function (resp. X π s ). Grey areas are estimates of the set ∆ μ (second row) and ∆μ (third row). Latent image correspond to nonparametric estimations of the Voronoi density estimator on W (first row), ∆ μ (second row) and ∆μ (third row). Titles include the observed number of points and the p-value of the quadrat test (denoted by QT).

Figure A. 9 :

 9 Figure A.9: First row present two simulations with on average 50,000 points of the sin model. Second row (resp. third row) present two subsampled versions X π s with the theoretical retaining probability function (resp. X π s ). Grey areas are estimates of the set ∆ μ (second row) and ∆μ (third row). Latent image correspond to nonparametric estimations of the Voronoi density estimator on W (first row), ∆ μ (second row) and ∆μ (third row). Titles include the observed number of points and the p-value of the quadrat test (denoted by QT).

  The graphical results together with the global combined p-value equal to 0.05% are unambiguous. The departure to the Poisson modeling is highly significant. In other words, there is strong evidence to reject the inhomogeneous Poisson model for the lightning strikes dataset. The two plots also show that K inh,t (τ ) > 2τ and K inh,s (r) > πr 2 for all values of r, τ considered and suggest that the X st is probably clustered a lot in time and space.

		Combined global test: p=0.005
		0.03				4
						3
	K t (τ)	0.02			K s (r)	2
		0.01
						1
		0.00				0
		0.000	0.002	0.004 τ	0.006	0.0	0.2	r	0.4	0.6
			Data	Central	Poisson		Data	Central	Poisson

Acknowledgements

This research is funded by Labex PERSYVAL-lab ANR-11-LABX-0025. The authors would like to take the opportunity to thank first Météorage and Méteo-France and especially Mickaël Zamo for providing us with the data and for scientific exchanges and second Jiri Dvorak for discussions and for providing us with some R code.

and Computing in Applied Probability, 16:411-431.

Appendix A. Homogeneous point process from a subsampled inhomogeneous point process In this section we present a simulation study in order to illustrate the subsampling strategy (b) presented in Section 2.4. This strategy intends to make an inhomogeneous point process a homogeneous one. We use a toy example in the planar case (thus no temporal component). Let λ s be the inhomogeneous spatial intensity of X s . We propose the following practical algorithm to implement strategy (b). Let ν π , set by the user, be the expected number of points of the subsampled version X π s .

(i) We estimate non parametrically the intensity λ s using the Voronoi-based nonparametric estimator introduced by Barr (iii) Let X π s be a subsampled version of X s obtained with retention probability π(x) = µ λs(x)

for any x ∈ ∆µ , thus defined on ∆µ ⊂ W .

The output of this algorithm is a realization of X π s whose intensity is close to µλ s (x)/ λs (x) ≈ µ for any x ∈ ∆µ . Thus X s is expected to be homogeneous. To illustrate this procedure, we consider two models for λ s in the planar case. We focus only on the Poisson case and test the homogeneity assumption of a pattern using the quadrat test procedure [START_REF] Baddeley | Spatial point patterns: methodology and applications with R[END_REF]