[R-sig-Geo] rainfall interpolation using ANN

Tue Dec 2 18:17:59 CET 2014

Peter,
This is a simplified version of your script untested, of course, given you
send you no data):

library(gstat)
library(raster)

dat71 <- shapefile("d:/FireDanger/IDW/71stations.shp")
sample.grid <- raster("rawsd_71Stations.tif")
#note that there are parameters you need to choose
m <- gstat(formula=VALUE1 ~1, data=dat71, nmax=10, set=list(idp = 0.5))
r_grid71 <- interpolate(sample.grid, m)

plot(r_grid71)

Robert

On Tue, Dec 2, 2014 at 8:45 AM, ping yang <pingyang.whu at gmail.com> wrote:

> Hi Nahm Lee and R-sig-geo,
>
> 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:
>
> library(gstat)
> library(maptools)
> library(raster)
> dat71 <- readShapePoints("d:/FireDanger/IDW/71stations.shp")
> bubble(dat71,"VALUE1")
> sample.grid <- raster("rawsd_71Stations.tif")
> Left <- bbox(sample.grid)[1]
> Right <- bbox(sample.grid)[3]
> Bottom <- bbox(sample.grid)[2]
> Top <- bbox(sample.grid)[4]
> cell.size <- res(sample.grid)
> grid.xy <- expand.grid(x = seq(Left, Right, cell.size[1]), y = seq(Top,
> Bottom, -cell.size[2]))
> coordinates(grid.xy) <- ~x + y
> gridded(grid.xy) = TRUE
> value71.idw <- krige(VALUE1 ~ 1, dat71, grid.xy)
> spplot(value71.idw["var1.pred"])
> r_grid71 <- raster(value71.idw["var1.pred"])
> plot(r_grid71)
>
> 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).
>
> Looking forward to hearing from you.
>
> Peter
>
> On Mon, Dec 1, 2014 at 12:31 PM, Nahm Lee <nlee at valleywater.org> wrote:
>
>>  This is my example with idw based on a web site, #
>> http://www.geo.ut.ee/aasa/LOOM02331/R_idw_interpolation.html.
>>
>>
>>
>> Nahm Lee
>>
>> nlee at valleywater.org
>>
>>
>>
>>
>>
>>
>>
>> # http://www.geo.ut.ee/aasa/LOOM02331/R_idw_interpolation.html
>>
>>
>>
>> library(ggplot2)
>>
>> library(gstat)
>>
>> library(sp)
>>
>> library(maptools)
>>
>> library(rgeos)
>>
>>
>>
>> setwd("C:/Users/Michael/Dropbox/R/SHP")
>>
>>
>>
>> AlertRain <- read.csv(file = "AlertRain.csv",header = TRUE)
>>
>>
>>
>> names(AlertRain)
>>
>>
>>
>> AlertRain_test <- AlertRain   # duplicate air temp. data file
>>
>> AlertRain$x <- AlertRain$longitude  # define x & y as longitude and
>> latitude
>>
>> AlertRain$y <- AlertRain$latitude
>>
>>
>>
>> # Set spatial coordinates to create a Spatial object:
>>
>> coordinates(AlertRain_test) = ~x + y
>>
>>
>>
>> plot(AlertRain_test)
>>
>>
>>
>> # Define the grid extent:
>>
>> #x.range <- as.numeric(c(21.76, 28.21))  # min/max longitude of the
>> interpolation area
>>
>> #y.range <- as.numeric(c(57.45, 59.72))  # min/max latitude of the
>> interpolation area
>>
>>
>>
>> x.range <-
>> c(min(round(AlertRain_test$x,2)-0.1),max(round(AlertRain_test$x,2)+0.1))
>>
>> y.range <-
>> c(min(round(AlertRain_test$y,2)-0.1),max(round(AlertRain_test$y,2)+0.1))
>>
>>
>>
>> # Create a data frame from all combinations of the supplied vectors or
>> factors.
>>
>> # See the description of the return value for precise details of the way
>> this is done.
>>
>> # Set spatial coordinates to create a Spatial object. Assign gridded
>> structure:
>>
>> grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1),
>> y = seq(from = y.range[1],
>>
>>     to = y.range[2], by = 0.1))  # expand points to grid
>>
>> coordinates(grd) <- ~x + y
>>
>> gridded(grd) <- TRUE
>>
>>
>>
>> #Plot the weather station locations and interpolation grid:
>>
>> plot(grd, cex = 1.5, col = "grey")
>>
>> points(AlertRain_test, pch = 1, col = "red", cex = 1)
>>
>>
>>
>> # interpolate surface and fix the output:
>>
>>
>>
>> idw <- idw(formula = Elevation ~ 1, locations = AlertRain_test,
>>
>>     newdata = grd)  # apply idw model for the data
>>
>> ## [inverse distance weighted interpolation]
>>
>>
>>
>> idw.output = as.data.frame(idw)  # output is defined as a data table
>>
>> names(idw.output)[1:3] <- c("long", "lat", "var1.pred")  # give names to
>> the modelled variables
>>
>>
>>
>> #Plot the results:
>>
>>
>>
>> ggplot() + geom_tile(data = idw.output, aes(x = long, y = lat, fill =
>> var1.pred)) + scale_fill_gradient(low = "cyan", high = "orange")+
>>
>>     geom_point(data = AlertRain, aes(x = longitude, y = latitude), shape
>> = 21,
>>
>>         colour = "red")+labs(title = "Elevation by rainfall station
>> elevation")
>>
>>
>>
>>
>>
>>
>>
>> *From:* R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] *On Behalf Of
>> *Mengxi Yang
>> *Sent:* Monday, December 01, 2014 1:32 AM
>> *To:* r-sig-geo at r-project.org
>> *Subject:* [R-sig-Geo] rainfall interpolation using ANN
>>
>>
>>
>> Dear list,
>>
>> I am going to do daily rainfall interpolation using ANN. but I am not
>> familiar with it. Can somebody give some advise or example?
>>
>>
>>
>> my data like:
>>
>>  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.
>>
>> [image: Inline image 2]
>>
>>
>>
>> can someone give me some advice? because I want to compare ANN with
>> Kriging and IDW.
>>
>> Thank you very much.
>>
>>
>>
>> Mengxi Yang
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20141202/71066e42/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 14586 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20141202/71066e42/attachment.png>