[R-sig-Geo] daily precipitation idw & kriging and beginner questions

[Paul Hiemstra]

Hi Ivan,

An interesting paper in this regard is:

Schuurmans, J.; Bierkens, M.; Pebesma, E. & Uijlenhoet, R. Automatic 
prediction of high-resolution daily rainfall fields for multiple 
extents: The potential of operational radar Journal of Hydrometeorology, 
2007, 8, 1204-1224

You are working at the KNMI, so you should be able to get a hold of the 
radar imagery for that day. I know that the author of the paper did all 
the work in R.

cheers,
Paul

Soenario, Ivan (KNMI) wrote:
> Hello experts,
>
>  
>
> Here's another novice to R. Hope you can help me out. Apologies in
> advance if questions have been asked before (in which case, could you
> tell me the thread).
>
>  
>
> My dataset is a table with daily Precipitation for the month January in
> 2000, measured at more than 300 weather stations. Each row is a weather
> station and there are 31 columns (or should I say attribute), each day
> it's own column.
>
>  
>
> I use a for-loop to go through the columns (the days) and each loop
> performs an idw and autokrige interpolation, plots images and saves them
> as jpg to hard disk.
>
>  
>
> Here's the script:
>
>  
>
> Pd200001 <- read.table(paste(wordir, "/Input/", sep=""), sep=("\t"),
> header = TRUE)
>
> coordinates(Pd200001) = ~RD_LOCATION_X + RD_LOCATION_Y
>
> nl.grd <- read.asciigrid(paste(wordir, "/Input/nl_clip_1km.asc",
> sep=("")), as.image=FALSE, plot.image=TRUE)
>
>  
>
> for (i in 1:31) {
>
> kol <- paste(names(Pd200001[i+2]), "~1", sep="")     # first two columns
> are id and stationnr 
>
>  
>
> # INVERSE DISTANCE
>
> Pd.idw <- krige(as.formula(kol), Pd200001, nl.grd)
>
> image(Pd.idw["var1.pred"], col = topo.colors(20))
>
> title(paste("Pd_idw", 20000100+i, sep = ""))
>
> dev.copy(jpeg, file = paste( wordir, "/img/Pd_idw_", 20000100+i, ".jpg",
> sep = ""), height = 17, width = 17, units = "cm", quality = 100,
> res=150, bg="white") ; dev.off()
>
> rm(Pd.idw)
>
>  
>
> # ORDINARY KRIGING
>
> Pd.kri <- autoKrige(as.formula(kol), Pd200001, nl.grd)
>
> plot(Pd.kri)
>
> dev.copy(jpeg, file = paste( wordir, "/img/Pd_kri_", 20000100+i, ".jpg",
> sep = ""), height = 17, width = 17, units = "cm", quality = 100,
> res=150, bg="white") ; dev.off()
>
> rm(Pd.idw)
>
> }
>
>  
>
> The loop works and I get images in files but I also get Warning
> Messages, both after interpolation and after plot().
>
>  
>
>   
>> Pd.idw <- krige(as.formula(kol), Pd200001, nl.grd)
>>     
>
> [inverse distance weighted interpolation]
>
> Warning message:
>
> In points2grid(points, tolerance, round, fuzz.tol) :
>
>   grid has empty column/rows in dimension 2
>
>  
>
>   
>> spplot(Pd.idw, "var1.pred")
>>     
>
> Warning message:
>
> In data.row.names(row.names, rowsi, i) :
>
>   some row.names duplicated: 2,3,4,5,6,7,8,9,10,
>
>  
>
> Question 1
>
> How "bad" are the Warning Messages, how should I prevent them
>
>  
>
> Question 2
>
> I use 2 different plotting techniques, which makes the images hard to
> compare. At this stage I do not need real map data for further analysis,
> just some jpegs to quickly compare the results. How can I plot the idw's
> and de krige results the same way? They are of different class.
>
>  
>
> Question 3
>
> automapPlot(plot_data, zcol, col.regions, ...)
>
> The ... is intriguing: "Arguments that are passed on to spplot." This is
> also the case for many R functions, which makes R very flexible, but for
> a beginner this seems endlessly layered. For instance, I would like to
> create my own colour ramp, in either image(), plot() or spplot(). But I
> can't figure it out. Is it also possible to reverse the standard ramps
> heat.colors, topo.colors etc?
>
>  
>
> Question 4
>
> After the jpg is written to HD, I remove the objects Pd.idw and Pd.kri.
> Need I do this, or is the object overwritten anyway, at the start of the
> next loop? Also, is there another method to deal with the results, for
> instance can I create an object Pd.idw.many which contains all of the 31
> interpolated maps? This would be convenient for further R processing.
>
>  
>
> Many thanks in advance!
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul