[prev in list] [next in list] [prev in thread] [next in thread]
List: r-sig-geo
Subject: Re: [R-sig-Geo] fitting variogram in regrssion kriging
From: Helen Chen <hc10024 () gmail ! com>
Date: 2014-09-26 9:27:07
Message-ID: CAKjCKvGj2GAxeebx-QnaeNRfvsskvGiAPEGBM=vq_frVW6UeAA () mail ! gmail ! com
[Download RAW message or body]
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@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@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@uni-muenster.de
> > <mailto:edzer.pebesma@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@r-project.org <mailto:R-sig-Geo@r-project.org>
> > >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >>
> > >
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > R-sig-Geo@r-project.org <mailto:R-sig-Geo@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@r-project.org <mailto:R-sig-Geo@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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic