[R-sig-Geo] fitting variogram in regrssion kriging
Helen Chen
hc10024 at gmail.com
Fri Sep 26 11:27:07 CEST 2014
Hi,
Thanks for your response. The mask_10m.asc is a DEM data, not a column of
precipitation.
The precipitation data is generated from 42 rain gauges data. The first
column is Y, the second column is X, and the third column is rainfall
amounts (mm).
The following is the script I adopted from the regression-Kriging guide
website:
library(sp)
library(lattice)
trellis.par.set(sp.theme()) # plots the final predictions using
blue-pink-yellow legend
precip <- read.csv("2001_stations_rainfall_0923.csv",header=T)
str(precip)
coordinates(precip)=~Lon+Lat # this makes depth a SpatialPointsDataFrame
str(precip)
elevation = read.asciigrid("mask_10m.asc") # reads ArcInfo Ascii raster map
eleok=na.omit(elevation)
elevation[!is.na(elevation),]
str(eleok)
spplot(eleok, scales=list(draw=T), sp.layout=list("sp.points", precip,
pch="+"))
#Plot the xy graph target versus predictor:
elev.ov = overlay(eleok, precip) # create grid-points overlay
str(elev.ov at data)
precip$mask_10m.asc =elev.ov$mask_10m.asc # copy the slope values
lm.precip <- lm(annaul_2001~mask_10m.asc, as.data.frame(precip))
summary(lm.precip)
plot(annaul_2001~mask_10m.asc, as.data.frame(precip))
abline(lm(annaul_2001~mask_10m.asc, as.data.frame(precip)))
#Fit the variogram model of the residuals:
library(gstat)
null.vgm <- vgm(var(precip$annaul_2001), "Sph",
sqrt(areaSpatialGrid(elevation))/4, nugget=0) # initial parameters
vgm_precip_r <- fit.variogram(variogram(annaul_2001~mask_10m.asc, precip),
model=null.vgm)
plot(variogram(annaul_total_2001.mm.~mask_10m.asc, precip), vgm_precip_r,
main="fitted by gstat")
I am stuck in the "vgm_precip_r <-
fit.variogram(variogram(annaul_2001~mask_10m.asc, precip), model=null.vgm)"
line
Thanks again for your time and help.
Helen
2014-09-26 2:17 GMT-07:00 Edzer Pebesma <edzer.pebesma at uni-muenster.de>:
> Is mask_10m.asc a component of (column in) precip, or is it a separate
> object? We need to see the full script that includes the steps that were
> taken to constructed precip in order to help.
>
> On 09/26/2014 09:35 AM, Helen Chen wrote:
> > Hi,
> >
> > Thanks for your response. I have used traceback() and found out that the
> > bug is in the mask_10m.asc.
> >
> > The following is the message:
> >
> > 12: stop("missing values in object")
> > 11: na.fail.default(list(annaul_2001 = c(945.134, 648.716, 559.562,
> > 655.574, 771.144, 805.434, 1276.604, 519.938, 1036.574, 763.016,
> > 757.428, 1062.482, 851.662, 720.852, 672.084, 855.218, 599.948,
> > 745.744, 685.292, 699.77, 775.97, 752.602, 573.786, 755.65, 628.904,
> > 849.63, 741.426, 676.402, 909.066, 695.452, 577.85, 707.644,
> > 762, 824.738, 919.48, 653.034, 781.304, 886.968, 828.802, 842.518,
> > 459.994, 624.84), mask_10m.asc = c(88.5786418914795,
> 45.8059759140015,
> > 15.0768859386444, 12.3127512931824, 71.3703689575195,
> 152.519027709961,
> > 1005.87916564941, 11.7274975776672, NA, 228.937419891357,
> > 181.906177520752,
> > 318.3203125, 502.518295288086, NA, 42.9922714233398, NA, NA,
> > 104.217571258545, 177.359172821045, 232.616687774658,
> 457.08145904541,
> > 398.026458740234, NA, 232.151866912842, 64.8636112213135,
> > 320.606918334961,
> > 154.726047515869, 240.452461242676, NA, 293.9072265625, NA,
> > 181.774570465088,
> > 222.223022460938, 192.939517974854, 583.289840698242,
> 782.529495239258,
> > 497.32640838623, 1249.03039550781, 347.13688659668, NA, NA, NA
> > )))
> >
> > I have tried to use na.omit to remove the NAs from mask_10m.asc, but it
> > seems not work in this case.
> >
> > Do you have any idea about how to get rid of NAs from ascii data? I have
> > attache the ascii file and precipitation file I used.
> >
> > Many thanks in advance.
> >
> > Best,
> >
> > Helen
> > mask_10m.asc
> > <
> https://docs.google.com/file/d/0BxOsJNA0fHDrZjBSZURXNTBNWWM/edit?usp=drive_web
> >
> >
> >
> > 2014-09-25 2:04 GMT-07:00 Edzer Pebesma <edzer.pebesma at uni-muenster.de
> > <mailto:edzer.pebesma at uni-muenster.de>>:
> >
> > We can only speculate, but mask_10m.asc may contain missing values,
> and
> > not refer to rainfall (directly).
> >
> > On 09/25/2014 09:57 AM, Tomislav Hengl wrote:
> > >
> > > Helen,
> > >
> > > Most likely you need to subset the missing values by using e.g.
> > > "precip[!is.na <http://is.na>(precip$annaul_2001),]", but it is
> > difficult to see if you
> > > maybe have a typo with variable names.
> > >
> > > Try using "traceback()" and check the amount of missing values by
> > using
> > > e.g. "summary(is.na <http://is.na>(precip$annaul_2001))" and/or
> > > "summary(complete.cases(precip))".
> > >
> > > HTH,
> > >
> > > T. (Tom) Hengl
> > > Researcher @ ISRIC - World Soil Information
> > > Url: http://www.wageningenur.nl/en/Persons/dr.-T-Tom-Hengl.htm
> > > Network: http://profiles.google.com/tom.hengl
> > > Publications:
> http://scholar.google.com/citations?user=2oYU7S8AAAAJ
> > >
> > >
> > > On 25-9-2014 7:57, Helen Chen wrote:
> > >> Hi,
> > >>
> > >> I am following the code provide from this link:
> > >>
> >
> http://spatial-analyst.net/wiki/index.php?title=Regression-kriging_guide)
> > >> with my own data.
> > >>
> > >> but I got stuck in fitting variogram, the error message is as
> follow:
> > >>
> > >>> vgm_precip_r <- fit.variogram(variogram(annaul_2001~mask_10m.asc,
> > >> precip), model=null.vgm)
> > >> Error in na.fail.default(list(annaul_2001 = c(945.134, 648.716,
> > >> 559.562, :
> > >> missing values in object
> > >>
> > >> But I have checked that there is no missing value in my rainfall
> > data.
> > >> Any
> > >> help would be very much appreciated.
> > >>
> > >> Helen
> > >>
> > >> [[alternative HTML version deleted]]
> > >>
> > >> _______________________________________________
> > >> R-sig-Geo mailing list
> > >> R-sig-Geo at r-project.org <mailto: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 <mailto:R-sig-Geo at r-project.org>
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> > --
> > Edzer Pebesma
> > Institute for Geoinformatics (ifgi), University of Münster
> > Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
> > 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795
> >
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
> 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795
>
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list