How to convert converted data

Question

How to convert converted data

I have daily rainfall from 61 calibration stations for 12 years in the catchment area (8000 km2).

The goal is to create 5 km and 25 km trellised daytime rain. Since the stations are not small, and not all stations have rain even in the rainy season, I choose to go with a climatological variogram.

The rainfall measurement for a typical day (irain) is shown below with a few missing values indicated by NA.

7.8  4.4 15.4 19.1  5.8    0   42  6.4   21    21     0     0     0  15.6 0     0    10     5   1.2     0  14.4    NA    25  13.2     0   9.2 2   6.6   7.8  13.2  15.4    NA     9     0  15.5     0  18.6     6 0   4.8  10.6     0    9     0  12.4    NA    12     0     3    14 10    10     0    68  21.8    18  14.8   5.4     7     0    NA

Since the daily rainfall is skewed, for the conversion I tested with cube root conversion and log conversion (log1p) for each day individually . However, both conversions are not suitable for all days, as I tested with the shapiro wilk test. So, I choose the normal conversion of the score counter (NCR), as suggested by Grimes and Pardo (2009) Geostatistical analysis of precipitation . and used by Prof. Ashton Shortridge Code

The code below is used to generate a climatological semivariogram for the monsoon season. Notice that I used days when more than 30% of the stations reported rain. I did not find links to this. Choose 30%, since I get about 65% of the days to pass the threshold.

lag = 3
bins.vario = seq(0, 75, lag)
nb.bins = length(bins.vario) - 1 

nb.classes = numeric(nb.bins)  
vario.emp = array(0,c(nb.bins,6))
variance.emp = array(NA,c(10000,nb.bins))
vario.emp = as.data.frame(vario.emp)
class(vario.emp) = c("gstatVariogram","data.frame")
names(vario.emp) = c("np","dist","gamma","dir.hor","dir.ver","id")


nRows = nrow(stn.rain.subset)
for (i in 1:nRows) 
{   
  irain = stn.rain.subset[i,]
  isMissing = is.na(irain)
  isZero = (irain == 0)
  irain = irain[!isMissing & !isZero]
  irain = as.numeric(irain)
  rain.mean[[i]] = mean(irain); rain.var[[i]] = var(irain);

  # Testing with log-transformation
  # irain.logtt = log1p(irain)
  # # qqPlot(irain.logtt,distribution = "norm")
  # res = shapiro.test(irain.logtt)
  # pval[[i]] = res$p.value

  if (var(irain)>0)
  {
    # Scaling of the rain on each day. rain in mm. After scaling this is unitless
    irain.scaled = irain/sqrt(var(irain))
    irain.nsc = nscore(irain)
    score  = irain.nsc$nscore   

    easting = lon.UTM[!isMissing & !isZero] # Removing the stn.locs with NA values
    northing = lat.UTM[!isMissing & !isZero]

    rain = data.frame(easting,northing,score)
    names(rain) = c("easting","northing","rain")
    coordinates(rain) = c("easting","northing")  
    proj4string(rain) = UTM

    v = variogram(rain~1, data = rain,boundaries = bins.vario)

    bin.nb = ceiling(v$dist/lag)
    nb.classes[bin.nb] = nb.classes[bin.nb]+1
    vario.emp[bin.nb,] = vario.emp[bin.nb,]+v

}

:

, .

"model" "psill" "range" "kappa" "ang1"  "ang2"  "ang3"  "anis1" "anis2"
"Nug"   0.446609415762287   0   0   0   0   0   1   1
"Sph"   0.533499909345213   51.7206027419321    0.5 0   0   0   1   1

Similary for kriging read every day non-zero, . , ?

rain = data.frame(easting,northing,score)
# Multiplying the nugget and sill with var(rain) for each day.
clim.vrmod$psill = clim.vrmod$psill * var(irain)
krige.ok = krige(rain[,3]~1,locations =~easting+northing,data = rain,newdata=output.grd,model = clim.vrmod) 
krige.ok$var1.pred.bt = (backtr(krige.ok$var1.pred,irain.nsc, tails='separate'))*sqrt(var(irain))
krige.ok$var1.se = (krige.ok$var1.var)

, , :