[R-sig-Geo] writing shapefiles using write.pointShape
modern82376 at mypacks.net
modern82376 at mypacks.net
Thu Mar 13 22:10:27 CET 2008
I've been running gwr on some survey data and trying to write output back out to point shapefiles. At first, I thought that the shapefile wasn't being written, but after calling the file object I found it in a local settings/temp folder with a random shapefile name. I thought I had specified all that, but I'm obviously not familiar enough with the syntax to get it right. Would someone mind assisting me on this? Code follows..
> library(spgwr)
> surveypts <- readShapePoints("C:/temp/NE_survey.shp", proj4string = CRS("+proj=aea +lat_0=96w +lat_1=29.5n +lat_2=45.5n"))
> survey.bw <- gwr.sel(GHDENSITY ~ pca1 + PCA2 + PCA3 + PCA4 + PCA5, data=surveypts)
> fitpts <- readShapePoints("C:/temp/fitpts.shp", proj4string = CRS("+proj=aea +lat_0=96w +lat_1=29.5n +lat_2=45.5n"))
> fitpix <- SpatialPixels(fitpts, proj4string = CRS("+proj=aea +lat_0=96w +lat_1=29.5n +lat_2=45.5n"))
> surv.local <- ggwr(GHDENSITY ~ pca1 + PCA2 + PCA3 + PCA4 + PCA5, data = surveypts, bandwidth = survey.bw, gweight = gwr.Gauss, adapt = NULL, fit.points = fitpix, longlat = FALSE)
> surv.local
Call:
ggwr(formula = GHDENSITY ~ pca1 + PCA2 + PCA3 + PCA4 + PCA5,
data = surveypts, bandwidth = survey.bw, gweight = gwr.Gauss,
adapt = NULL, fit.points = fitpix, longlat = FALSE)
Kernel function: gwr.Gauss
Fixed bandwidth: 22738.78
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. -4.176e+03 -9.067e+01 -3.565e+00 8.029e+01 4.402e+03 -73.5570
pca1 -2.567e+00 -2.289e-02 -1.725e-03 1.674e-02 3.135e-01 -0.0077
PCA2 -7.874e+01 -6.543e-01 -1.934e-01 2.298e-01 2.016e+01 -0.2667
PCA3 -8.844e+01 -1.236e+00 -1.109e-01 6.008e-01 1.469e+02 0.0202
PCA4 -1.272e+01 -3.509e-01 3.005e-01 1.178e+00 3.386e+02 0.8161
PCA5 -1.368e+02 -1.134e+00 -8.233e-03 9.222e-01 4.384e+01 0.6937
> names(surv.local$SDF)
[1] "sum.w" "X.Intercept." "pca1" "PCA2" "PCA3"
[6] "PCA4" "PCA5" "dispersion" "response_resids"
> coords <- coordinates(fitpts)
> pca1 <- cbind(coords, surv.local$SDF$pca1)
> summary(pca1)
coords.x1 coords.x2 V3
Min. :-670120 Min. :287533 Min. :-2.566983
1st Qu.:-544120 1st Qu.:375533 1st Qu.:-0.022891
Median :-417620 Median :464533 Median :-0.001725
Mean :-417620 Mean :464533 Mean :-0.047442
3rd Qu.:-291120 3rd Qu.:553533 3rd Qu.: 0.016741
Max. :-165120 Max. :641533 Max. : 0.313544
>class(pca1)
[1] "matrix"
> dim(pca1)
[1] 179630 3
> all(dim(dpca1 <- as.data.frame(pca1)) == c(179630,3))
[1] TRUE
> summary(dpca1)
coords.x1 coords.x2 V3
Min. :-670120 Min. :287533 Min. :-2.566983
1st Qu.:-544120 1st Qu.:375533 1st Qu.:-0.022891
Median :-417620 Median :464533 Median :-0.001725
Mean :-417620 Mean :464533 Mean :-0.047442
3rd Qu.:-291120 3rd Qu.:553533 3rd Qu.: 0.016741
Max. :-165120 Max. :641533 Max. : 0.313544
> class(dpca1)
[1] "data.frame"
> B1pca1 <- tempfile("C:\temp")
> write.pointShape(coords, df=dpca1, B1pca1, factor2char=TRUE)
(Where's the shapefile??)
>B1pca1
[1] "C:\\DOCUME~1\\gcook\\LOCALS~1\\Temp\\Rtmpijon4d\\C:\temp1eb26e9"
After calling B1pca1, I find a randomly named shapefile (with .dbf and .shx)extensions in the .../Local Settings/Temp/ folder. I'm trying to name the file "B1pca1.shp" and I want it to go in the C:/Temp folder (for now). I've tried a few permutations of the last two commands, but I'm just not getting it right.
Thanks so much :)
Gericke
More information about the R-sig-Geo
mailing list