Re: [R-sig-Geo] Clipping a Map
Hi Ravi, After performing 'gIntersects' try merging this result to your dataframe, then select the TRUE values? Jesse cz_zip - gIntersects(camapzip_temp,camap_base, byid=TRUE) camap_base@data-cbind(camap_base@data, czip) camap2-camap_base[(camap_base@data[,/???/] %in% c(TRUE)),] # Insert your column number for czip in place of ??? plot(camap2) head(camap2@data) - Jesse D Berman, PhD Yale University School of Forestry and Environmental Studies Post-Doc Fellow -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Clipping-a-Map-tp7585156p7585162.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot class intervals
Milan, Try using the function 'cuts' instead of at=c(). Jesse library(sp) data(meuse) coordinates(meuse)=~x+y spplot(meuse, cadmium, do.log = TRUE, key.space=list(x=0.2,y=0.9,corner=c(0,1)), cuts=c(1, 2, 3, 4, 10, 20), # My custom cuts scales=list(draw=T)) - Jesse D Berman, PhD Yale University School of Forestry and Environmental Studies Post-Doc Fellow -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/spplot-class-intervals-tp7584934p7584935.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] overlay a ppp object and a SpatialPolygon
Hi Laetitia, Try turning your .ppp object into a spatial grid data frame, then convert it to a raster. The 'extract' function should then work if your projections are correct. Jesse your.ppp.object-as.SpatialGridDataFrame.im(your.ppp.object) rast.ppp.object- raster(your.ppp.object) #Convert to a raster grid output- extract(rast.ppp.object, polys) - Jesse D Berman, PhD Yale University School of Forestry and Environmental Studies Post-Doc Fellow -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/overlay-a-ppp-object-and-a-SpatialPolygon-tp7584755p7584775.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How to do Zonal Statistics after Kriging in R given Shapefile of Polygon
James, it seems like your projections for 'town' and 'muni.sp' are different with one being lat/long and the other being a UTM Zonal. Try setting them the same and see if that solves your problem. Jesse - Jesse D Berman, PhD Yale University School of Forestry and Environmental Studies Post-Doc Fellow -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/How-to-do-Zonal-Statistics-after-Kriging-in-R-given-Shapefile-of-Polygon-tp7584330p7584595.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Singular Matrix Error
Hi All, I want to perform ordinary kriging on a series of air pollution values, but am getting the error / IL.Jul7.ok-krige(Arithmetic.Mean~1, IL.Jul7, model=ILvgm.jul7, newdata=IL2.map) [using ordinary kriging] chfactor.c, line 131: singular matrix in function LDLfactor() Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim, : LDLfactor/ I know that this error is typically caused by either duplicate points or a variogram without a nugget effect. However, I have neither (see code below). I also have set the same projection for points and prediction grid, plus ran an IDW successfully /[CRS(+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs)] / / ILvgm.jul7 model psill range 1 Nug 1.537923 0.000 2 Sph 53.045879 266.323 zerodist(IL.Jul7) [,1] [,2] summary(IL.Jul7@data$Arithmetic.Mean) Min. 1st Qu. MedianMean 3rd Qu.Max. 13.50 17.60 21.95 22.82 27.83 32.10 / Any thoughts as to what might be causing this problem? Sorry for not providing a reproducible example, as I can't reproduce it and have run similar code without problems on other data. Thanks for the suggestions. Jesse - Jesse D Berman, PhD Johns Hopkins Bloomberg School of Public Health Department of Environmental Health Sciences Post-Doc -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Singular-Matrix-Error-tp7584203.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] predict function of regressors in the raster package
Hi Hein, I'm not sure if this will help, but one thing to check is that your prediction grid has covariate data for each of the 40,000 cells. If a large number of cells have 'NA' as data values, then sometimes the prediction will not work. Offhandedly, it strikes me that ycoord may be limited to only your t0 locations. Jesse -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/predict-function-of-regressors-in-the-raster-package-tp7583201p7583205.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Maptools package
Hi all, I've loaded the new version of R 3.0.0 and it seems the package 'maptools' is not available on the mirrors. It's not urgent, but I'm curious if anyone knows when it might become available? Thanks. Regards, Jesse -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Maptools-package-tp7583198.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Large Prediction Variances with gstat
Hi All,
First time post, so please excuse any omissions/confusion. I am performing
a series of prediction models using gstat and discovered that prediction
variance of spatially dependent data with OLS models was larger than those
of kriging models. This is counter-intuitive to the assumption that
treating spatially dependent data as IID will result in artificially
shrunken prediction variances. Anyhow, to better understand how gstat
treats OLS predictions, I reproduced an OLS prediction with the base package
('predict' and 'predict.lm') and found that while I got identical beta's, I
got substantially higher variances with gstat.
Can anyone shed some light as to why gstat might be giving these larger
prediction variances when performing an OLS model? (see reproducible
example below)
Regards and thanks for the help,
Jesse
library(gstat)
data(meuse)
coordinates(meuse) = ~x+y
meuse.ns-as.data.frame(meuse) #non-spatial Meuse data
#Data sets for modeling and prediction; spatial and non-spatial
dat.mod-meuse[1:100,]
dat.pred-meuse[101:155,]
dat.ns.mod-meuse.ns[1:100,]
dat.ns.pred-meuse.ns[101:155,]
#Linear Model Prediction with base package
lm.zn-lm(log(zinc)~x+y+elev, data=dat.mod)
lm.zn.pred-predict(lm.zn, dat.pred, se.fit=TRUE)
pred.variance-(lm.zn.pred$se.fit)^2
#Linear Model Prediction with gstat
gstat.lm.zn.pred-krige(log(zinc)~x+y+elev, dat.mod, newdata=dat.pred)
#Compare Results
summary(gstat.lm.zn.pred@data) #gstat results
summary(lm.zn.pred$fit)
summary(pred.variance) #base model prediction variance
--
View this message in context:
http://r-sig-geo.2731867.n2.nabble.com/Large-Prediction-Variances-with-gstat-tp7582882.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] crop and overlay question
Hi Erin, I'm sure someone else will have a more elegant solution, but if you need a quick fix try this (adapted from your posted code) Jesse tmp - map(state,Texas,fill=TRUE,plot=FALSE) texas.boundary - map2SpatialPolygons(tmp, IDs=tmp$names, proj4string=CRS(+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0)) texas1.boundary - spTransform(texas.boundary,CRS(+proj=aea +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=150 +y_0=600 +ellps=GRS80 +datum=NAD83 +units=m +no_defs) ) grd - SpatialPoints(makegrid(texas1.boundary, n=300)) # New steps Overlay-overlay(grd, texas1.boundary) #To find out which points fall within the TX boundary grd@coords-cbind(grd@coords,Overlay) #To mark those points outside TX with an NA tex.grid-na.exclude(as.data.frame(grd)) # Transfers to non-spatial dataframe to remove NA points coordinates(tex.grid)=~x1+x2 # Now turn back into spatial points dataframe proj4string(tex.grid)=(+proj=aea +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=150 +y_0=600 +ellps=GRS80 +datum=NAD83 +units=m +no_defs) plot(tex.grid) -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/crop-and-overlay-question-tp7582894p7582896.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
