[R-sig-Geo] R: individuals distribution points in paired polygons
Alice C. Hughes
dr_achughes at hotmail.co.uk
Sun Sep 7 10:31:36 CEST 2014
I have distribution points (in excel or spatial point form) for a set ofindividuals and a file of polygons-with matching IDs: so each individual hasa set of points and a polygon with the same IDs: but there are 100's to1000's of individuals. What I would like to do is to go through and for each individual find thetotal number of points and the number which fall within the appropriatepolygon, then have the output in spreadsheet form with a column for ID;Total number of points; points in polygon. Each individual has a polygon with the same name, but every code I've seenso far can't differentiate, or scroll through multiple point files withunique IDs
Alice C. HughesAssociate ProfessorCentre for Integrative Conservation,Xishuangbanna Tropical Botanical Garden,Chinese Academy of SciencesMenglun, Mengla, Yunnan 666303, P.R. ChinaPh: 15198676559achughes at xtbg.ac.cn
> Date: Fri, 16 May 2014 14:32:22 +0200
> From: rubenfcasal at gmail.com
> To: mfiori at arpa.sardegna.it
> CC: r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] R: R: Regression Kriging cross validation
>
> Hi Michele,
>
> Kriging methods assume that the variogram is known, although it is not
> in practice. With your procedure, the whole prediction method is
> reproduced (adding the extra-variability due to the estimation of the
> variogram), so a better approximation to the variability of the
> prediction error should be (in principle) obtained.
>
> Excuse the rush...
>
> Ruben Fernandez-Casal
>
>
> El 16/05/2014 9:50, Michele Fiori escribió:
> > Dear All,
> > Finally I managed to develop this ad hoc procedure for the regression kriging cross validation; I tried to test the kriging part (the regression part has been already tested successfully) by comparison with the krige.cv function and I noticed that there are small discrepancies. I verified that they are due to the fact that in my case the variogram is recalculated for each step, while the function runs using the same initial variogram obtained on the complete set of data.
> > What do you think about?
> > Michele
> >
> > ARPAS - Environmental Protection Agency of Sardinia MeteoClimatic
> > Department - Meteorological Service
> >
> > Viale Porto Torres 119 - 07100 Sassari, Italy Tel + 39 079 258617 Fax
> > + 39 079 262681 www.sardegnaambiente.it/arpas
> >
> > #################### Cross validation REGRESSION KRIGING (leave on out)
> >
> > PP03.lm <- lm(reg, prec2)
> > Nstazioni <- length(prec2$NOME_STAZ)
> > for (i in 1:Nstazioni)
> > {
> > ### regression step for point “i”
> >
> > prec2.i <- prec2[-i,]
> > md.i <- lm(reg, data=prec2.i)
> > ypredlm.i <- predict(md.i,prec2[i,])
> > err[i] <- (ypredlm.i - prec2$PP03[i])^2
> > res[i] <- (ypredlm.i - prec2$PP03[i])
> >
> > ### kriging step for residual point “i”
> >
> > PP03i.rsvar <- variogram(prec2.i$residuals~1, prec2.i)
> > PP03i.ivgm <- vgm(nugget=0, model="Sph", range=sqrt(diff(prec2.i at bbox["x",])^2 + diff(prec2.i at bbox["y",])^2)/4, psill=var(prec2.i$residuals))
> > PP03i.rvgm <- fit.variogram(PP03i.rsvar, model=PP03i.ivgm)
> > ypredok.i<- krige(md.i$residuals~1, prec2.i, prec2[i,], PP03i.rvgm)
> > ypredok.i$var1.pred
> >
> > ### Regression + kriging predicted value for point “i”
> >
> > ypred.i <- ypredlm.i + ypredok.i$var1.pred
> >
> > err[i] <- (ypred.i - prec2$PP03[i])^2
> > ypred[i] <- ypred.i
> > }
> > err <-c(err)
> > ypred <-c(ypred)
> >
> > ### mean squared error
> >
> > MSE <- sum(err)/(Nstazioni)
> > MSE
> >
> > ### root mean square error
> >
> > RMSE.rk <- MSE^0.5
> > RMSE.rk
> >
> > Da: Moshood Agba Bakare [mailto:bakare at ualberta.ca]
> > Inviato: giovedì 15 maggio 2014 17:39
> > A: Michele Fiori
> > Cc: rubenfcasal; r-sig-geo at r-project.org
> > Oggetto: Re: [R-sig-Geo] R: Regression Kriging cross validation
> >
> > Hi Michele,
> > I have similar problem. I used ordinary kriging and inverse distance weighting method (IDW) to generate set of interpolated values from the same interpolation grid. I don't understand how cross validation can be done to come up with diagnostic statistics such as mse, rmse to use as basis for identifying the best interpolation method. I used krige.cv but I encountered error message. Please any advice on what to do please?
> >
> > ## Create grid for the interpolation through ordinary kriging and idw
> >
> > grid <- expand.grid(easting = seq(from = 299678, to = 301299, by=10),
> > northing=seq(from = 5737278, to = 5738129, by=10))
> >
> >
> > ## convert the grid to SpatialPixel class to indicate gridded spatial data
> >
> > coordinates(grid)<-~easting + northing
> > proj4string(grid)<-CRS("+proj=utm +zone=12 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")
> > gridded(grid)<- TRUE
> >
> > #### Ordinary kriging
> >
> > prok <- krige(id="yield",yield ~ 1, canmod.sp, newdata = grid, model=exp.mod,nmax=20,maxdist=33.0)
> >
> > ## Inverse Distance Weighting (IDW) Interpolation method
> >
> > yield.idw = idw(yield~1, canmod.sp, grid,nmax=20,maxdist=33.0,idp=1)
> >
> >
> > Thanks
> > Moshood
> >
> >
> > On Thu, May 15, 2014 at 9:23 AM, Michele Fiori <mfiori at arpa.sardegna.it> wrote:
> > Thank you for your kind reply
> > therefore as I have used the Osl method for regression, my result will never
> > match the universal kriging; However, in order to validate my method, I'm
> > trying to implement in the script a calculation loop witch runs n times (the
> > number of stations) regression + kriging without one station at a time.
> > Thank you again
> > Michele
> >
> >
> > -----Messaggio originale-----
> > Da: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org]
> > Per conto di rubenfcasal
> > Inviato: martedì 29 aprile 2014 19:49
> > A: r-sig-geo at r-project.org
> > Oggetto: Re: [R-sig-Geo] Regression Kriging cross validation
> >
> > Hello Michele,
> >
> > Universal kriging is equivalent to Linear Regression (with the
> > generalized-least-squaresestimator) + Simple Kriging of residuals (e.g.
> > Cressie, 1993, section 3.4.5). The differences you observe are probably due
> > to the use of ordinary least squares. If you use (leave-one-out)
> > cross-validation with krige.cv (considering the UK model), the trend is also
> > re-estimated at each prediction location. From my point of view, this would
> > be the recommended way to proceed.
> > As far as I know, there are no available implementations of the
> > procedure you are suggesting.
> >
> > Best regards,
> > Rubén.
> >
> >
> > El 29/04/2014 13:33, Michele Fiori escribió:
> >> Hi everyone,
> >> I am working on rainfall interpolation using regression kriging method
> >> and I need suggestions on how I can carry out a cross validation
> >> (leave-one-out) for this elaboration. At first I tried to apply
> >> directly Krige.cv, similarly to UK method (example for october:
> >> PP10uk.cv <- krige.cv(reg, prec2, PP10.vgm)), but unfortunately when I
> >> applied Universal Kriging on the same data, I realized that UK map was a
> > little different from RK map.
> >> So my question is: How could I manage universal kriging in order to
> >> make it equivalent to regression kriging and use the above
> >> cross-validation, or is there another different method to apply cross
> >> validation (leave-one-out) on Regression Kriging interpolation?
> >> Below my code:
> >> Many thanks
> >>
> >> Michele Fiori
> >>
> >> ARPAS - Environmental Protection Agency of Sardinia MeteoClimatic
> >> Department - Meteorological Service
> >>
> >> Viale Porto Torres 119 - 07100 Sassari, Italy Tel + 39 079 258617 Fax
> >> + 39 079 262681 www.sardegnaambiente.it/arpas
> >>
> >> #### Creating SpatialPixelDataFrame ("dem" - 250x250 m grid)
> >> ....
> >> #### Loading Precipitation data
> >> prec2 <- read.table("prec2.txt", sep="\t", header =TRUE)
> >> coordinates(prec2) <- c("x", "y")
> >> proj4string(prec2) <- CRS("+init=epsg:32632")
> >> #### Linear regression Model
> >> mod.gen <- lm(PP10 ~ QUOTA_MARE + UTM_EST + UTM_NORD + DIST_MARE,
> >> prec2)
> >> step1 <- stepAIC(mod.gen, direction="both")
> >> reg <- formula(step1)
> >> PP10.lm <- lm(reg, prec2)
> >> summary(PP10.lm)
> >> prec2$residuals <- residuals(PP10.lm)
> >> dem$predlm <- predict(PP10.lm, dem)
> >> #### Variogram of residuals
> >> PP10.vgm <- vgm(nugget=51.46, model="Sph", range=38038.89,
> >> psill=86.44)
> >> #### Ordinary Kriging of residuals
> >> PP10.okr <- krige(PP10.lm$residuals ~ 1, prec2, dem, PP10.vgm,
> >> maxdist=Inf)
> >> dem$varokr <- PP10.okr$var1.pred
> >> #### Regression Kriging (Linear Regression + Ordinary Kriging of
> >> residuals)
> >> dem$vark <- dem$predlm + dem$varokr
> >> #### Universal kriging
> >> PP10.uk <- krige(reg, prec2, dem, PP10.vgm, maxdist=Inf)
> >> dem$varuk <- PP10.uk$var1.pred
> >> dem$difference <- dem$vark - dem$varuk
> >> spplot(dem[c("difference")], col.regions=terrain.colors(25),
> >> contour=FALSE, cuts = 15)
> >>
> >> _______________________________________________
> >> 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
> >
> > _______________________________________________
> > 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
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list