<div dir="ltr">Peter, <div>This is a simplified version of your script untested, of course, given you send you no data):</div><div><br></div><div><div>library(gstat)<br></div><div>library(raster)</div><div><br></div><div>dat71 <- shapefile("d:/FireDanger/IDW/71stations.shp")</div><div>sample.grid <- raster("rawsd_71Stations.tif")</div><div>#note that there are parameters you need to choose<br></div><div>m <- gstat(formula=VALUE1 ~1, data=dat71, nmax=10, set=list(idp = 0.5)) </div><div>r_grid71 <- interpolate(sample.grid, m)</div><div><br></div><div>plot(r_grid71)</div><div><br></div><div>Robert</div><div><br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Dec 2, 2014 at 8:45 AM, ping yang <span dir="ltr"><<a href="mailto:pingyang.whu@gmail.com" target="_blank">pingyang.whu@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-family:georgia,serif;font-size:small">Hi Nahm Lee and R-sig-geo,<br><br>I am quite interested in the IDW interpolation comparing to other methods such as Krige, I have a code like this which I think also did the IDW:<br><br>library(gstat)<br>library(maptools)<br>library(raster)<br>dat71 <- readShapePoints("d:/FireDanger/IDW/71stations.shp")<br>bubble(dat71,"VALUE1")<br>sample.grid <- raster("rawsd_71Stations.tif")<br>Left <- bbox(sample.grid)[1]<br>Right <- bbox(sample.grid)[3]<br>Bottom <- bbox(sample.grid)[2]<br>Top <- bbox(sample.grid)[4]<br>cell.size <- res(sample.grid)<br>grid.xy <- expand.grid(x = seq(Left, Right, cell.size[1]), y = seq(Top, Bottom, -cell.size[2]))<br>coordinates(grid.xy) <- ~x + y<br>gridded(grid.xy) = TRUE<br>value71.idw <- krige(VALUE1 ~ 1, dat71, grid.xy)<br>spplot(value71.idw["var1.pred"])<br>r_grid71 <- raster(value71.idw["var1.pred"])<br>plot(r_grid71)<br><br>However, I am not sure it is correct. Have someone using this method to do the IDW? and One more question, How can I complete a Krige application based on a shapefile (there are other examples using the spreedsheet).<br><br>Looking forward to hearing from you.<br><br>Peter<br></div></div><div class="gmail_extra"><br><div class="gmail_quote"><div><div class="h5">On Mon, Dec 1, 2014 at 12:31 PM, Nahm Lee <span dir="ltr"><<a href="mailto:nlee@valleywater.org" target="_blank">nlee@valleywater.org</a>></span> wrote:<br></div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5">
<div link="blue" vlink="purple" lang="EN-US">
<div>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">This is my example with idw based on a web site, #
<a href="http://www.geo.ut.ee/aasa/LOOM02331/R_idw_interpolation.html" target="_blank">http://www.geo.ut.ee/aasa/LOOM02331/R_idw_interpolation.html</a>.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">Nahm Lee<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><a href="mailto:nlee@valleywater.org" target="_blank">nlee@valleywater.org</a><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># <a href="http://www.geo.ut.ee/aasa/LOOM02331/R_idw_interpolation.html" target="_blank">http://www.geo.ut.ee/aasa/LOOM02331/R_idw_interpolation.html</a><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">library(ggplot2)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">library(gstat)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">library(sp)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">library(maptools)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">library(rgeos)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">setwd("C:/Users/Michael/Dropbox/R/SHP")<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">AlertRain <- read.csv(file = "AlertRain.csv",header = TRUE)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">names(AlertRain)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">AlertRain_test <- AlertRain # duplicate air temp. data file<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">AlertRain$x <- AlertRain$longitude # define x & y as longitude and latitude<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">AlertRain$y <- AlertRain$latitude<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># Set spatial coordinates to create a Spatial object:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">coordinates(AlertRain_test) = ~x + y<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">plot(AlertRain_test)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># Define the grid extent:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">#x.range <- as.numeric(c(21.76, 28.21)) # min/max longitude of the interpolation area<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">#y.range <- as.numeric(c(57.45, 59.72)) # min/max latitude of the interpolation area<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">x.range <- c(min(round(AlertRain_test$x,2)-0.1),max(round(AlertRain_test$x,2)+0.1))<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">y.range <- c(min(round(AlertRain_test$y,2)-0.1),max(round(AlertRain_test$y,2)+0.1))<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># Create a data frame from all combinations of the supplied vectors or factors.
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># See the description of the return value for precise details of the way this is done.
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># Set spatial coordinates to create a Spatial object. Assign gridded structure:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1],
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"> to = y.range[2], by = 0.1)) # expand points to grid<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">coordinates(grd) <- ~x + y<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">gridded(grd) <- TRUE<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">#Plot the weather station locations and interpolation grid:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">plot(grd, cex = 1.5, col = "grey")<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">points(AlertRain_test, pch = 1, col = "red", cex = 1)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"># interpolate surface and fix the output:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">idw <- idw(formula = Elevation ~ 1, locations = AlertRain_test,
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"> newdata = grd) # apply idw model for the data<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">## [inverse distance weighted interpolation]<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">idw.output = as.data.frame(idw) # output is defined as a data table<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">names(idw.output)[1:3] <- c("long", "lat", "var1.pred") # give names to the modelled variables<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">#Plot the results:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d">ggplot() + geom_tile(data = idw.output, aes(x = long, y = lat, fill = var1.pred)) + scale_fill_gradient(low = "cyan", high = "orange")+<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"> geom_point(data = AlertRain, aes(x = longitude, y = latitude), shape = 21,
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"> colour = "red")+labs(title = "Elevation by rainfall station elevation")<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1f497d"><u></u> <u></u></span></p>
<div style="border:none;border-top:solid #b5c4df 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif"">From:</span></b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif""> R-sig-Geo [mailto:<a href="mailto:r-sig-geo-bounces@r-project.org" target="_blank">r-sig-geo-bounces@r-project.org</a>]
<b>On Behalf Of </b>Mengxi Yang<br>
<b>Sent:</b> Monday, December 01, 2014 1:32 AM<br>
<b>To:</b> <a href="mailto:r-sig-geo@r-project.org" target="_blank">r-sig-geo@r-project.org</a><br>
<b>Subject:</b> [R-sig-Geo] rainfall interpolation using ANN<u></u><u></u></span></p>
</div><div><div>
<p class="MsoNormal"><u></u> <u></u></p>
<div>
<p class="MsoNormal">Dear list,<u></u><u></u></p>
<div>
<p class="MsoNormal">I am going to do daily rainfall interpolation using ANN. but I am not familiar with it. Can somebody give some advise or example?<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">my data like:<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"> I have 325 stations in different location.each station have unique code and coordinates. And I want to interpolate in 1km resolution. I have the mask map of NL.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><img src="cid:image001.png@01D00D52.07333460" alt="Inline image 2" border="0" height="289" width="425"><u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">can someone give me some advice? because I want to compare ANN with Kriging and IDW.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal">Thank you very much.<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal">Mengxi Yang<u></u><u></u></p>
</div>
</div>
</div></div></div>
</div>
<br></div></div>_______________________________________________<br>
R-sig-Geo mailing list<br>
<a href="mailto:R-sig-Geo@r-project.org" target="_blank">R-sig-Geo@r-project.org</a><br>
<a href="https://stat.ethz.ch/mailman/listinfo/r-sig-geo" target="_blank">https://stat.ethz.ch/mailman/listinfo/r-sig-geo</a><br>
<br></blockquote></div><br></div>
<br>_______________________________________________<br>
R-sig-Geo mailing list<br>
<a href="mailto:R-sig-Geo@r-project.org">R-sig-Geo@r-project.org</a><br>
<a href="https://stat.ethz.ch/mailman/listinfo/r-sig-geo" target="_blank">https://stat.ethz.ch/mailman/listinfo/r-sig-geo</a><br>
<br></blockquote></div><br></div>