non-stationary precipitation model

[TheCorinna1994]

Dear all,

I am trying to implement a non-stationary precipitation model, where both the mean and dependence are allowed to vary with elevation. Unfortunately I dont know how to implement this. 

The only "dummy code" I have is from https://becarioprecario.bitbucket.io/spde-gitbook/ch-nonstationarity.html where they implement non-stationarity through the range by defining it as a function of the first coordinate of the location. 

The other code I have is from the R-inla tutorial. There they define a model where the local precision depends on the first coordinate.

Is there any other tutorial with code where the mean and dependence are allowed to vary with elevation? Why do they relate non-stationarity to the first coordinate?  Is there any preferred /right way to implement non-stationarity?

Your help will be highly appreciated. Thanks in advance already.

All the best,

Corinna

[Finn Lindgren]

Hi,

have you had a look at the Ingebrigtsen et al paper from ca 2014, that
discusses nonstationary precipitation modelling?

For the code examples, having it depend on the first coordinate is
just a way to get a basic covariate that is spatially varying (and
also not a completely unrealistic thing to have in many regions, with
winds blowing predominantly from west to east, making the "easting"
information relevant.

Finn

> --
> You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/b46ea051-553b-4f9a-83cd-49416c346ff1n%40googlegroups.com.



--
Finn Lindgren
email: finn.l...@gmail.com

[TheCorinna1994]

Dear Finn,

thank you for your answer. It helped a lot.

Yes I had a look at the Ingebrigtsen paper before, but unfortunately I do not really know how to implement 

log(tau(u))= theta_1^(tau)+ h(u)*theta_h^(tau) and log(kappa(u))=  theta_1^(kappa)+ h(u)*theta_h^(kappa) correctly.

I tried it though, where n=48 are my measurement stations, stored in est.data. 

theta_1.tau<-rnorm(n,mean= 4, sd=0.1)
theta_1.kappa<-rnorm(n, -4, 0.1)

theta_elev.tau<-rnorm(n,0,1)
theta_elev.kappa<-rnorm(n,0,1)

spde.non.stat.elev<-inla.spde2.matern(mesh,      # 1243 mesh vertices
                                      B.tau=cbind(0, theta_1.tau,0, est.data$Elevation*theta_elev.tau),
                                      B.kappa=cbind(0,0, theta_1.kappa, est.data$Elevation*theta_elev.kappa),
                                      theta.prior.mean=rep(0,3),
                                      theta.prior.prec=rep(1,3))

However I get the following error " B matrix has 48 rows but should have 1 or 1243 ".

Do you know how to fix that problem?

Thanks in advance already.

All the best,

Corinna

[Finn Lindgren]

Hi,
you need to have the covariate need for the non-stationarity
everywhere on the mesh, not just where you have observations of your
response variable.
Finn

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/802d647d-70b3-4e39-9f97-87b8725fb623n%40googlegroups.com.

[TheCorinna1994]

Hi Finn, 

thanks for your answer.

In my region of interest I have limited information, only 48 measurement stations with Precipitation, Elevation and Wind measurements.

Is there any way how to get Elevation Data for the mesh vertices or am I not able to do the non-stationary modelling then?

Best,

corinna

[TheCorinna1994]

P.S: not having covariate information on my mesh vertices in general, is a huge limitation for various modelling approaches expect for stationary modelling right? How did you get covariate information on mesh vertices?

Best,

Corinna

[Finn Lindgren]

For elevation, it depends on what spatial resolution is needed, but the publicly available ETOPO databases has served us well in the past. It’s a global dataset available in different resolutions. Haven’t looked at it in a while, so unsure of the options.

Finn

On 22 Mar 2023, at 14:35, TheCorinna1994 <corinn...@hotmail.com> wrote:

Hi Finn, 

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/90e53f66-2e9e-4ee7-95b5-a497971a4e36n%40googlegroups.com.

[Finn Lindgren]

I haven’t checked our old paper, but we either used ETOPO or had elevation info from some Norwegian source.

Finn

On 22 Mar 2023, at 14:43, TheCorinna1994 <corinn...@hotmail.com> wrote:

P.S: not having covariate information on my mesh vertices in general, is a huge limitation for various modelling approaches expect for stationary modelling right? How did you get covariate information on mesh vertices?

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/01422554-a2fe-47ce-b963-ad3cff27bfd3n%40googlegroups.com.

[TheCorinna1994]

Ok Ok thanks.

One last question: I have my mesh and now I just imported a file containing the elevation throughout entire Austria via the raster package.

How do I relate the mesh vertices to the altitude data frame?

[manjagoc]

Hello to everyone,

I have the same doubt. Can someone can give an example or refer me to a previous answer in this group?

[Finn Lindgren]

Hi,

I had missed the previous question.

When the covariate is on a finer scale than the mesh, I tend to do this:

A <- inla.spde.make.A(mesh, loc=covariategridlocations)

Covariateonmesh<-(Covariatevalues %*% A) / colSums(A)

This evaluates locally weighted averages of the covariates, using the mesh basis functions as weights.

Finn

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/3cd90a89-e5b6-41ec-91b4-8f8da3dd2b41n%40googlegroups.com.

[manjagoc]

Thank you very much Professor Lindgren.