[R-sig-Geo] writing shapefiles using write.pointShape

Sun Mar 16 11:42:17 CET 2008

On Thu, 13 Mar 2008, modern82376 at mypacks.net wrote:

> I'm also having issues with getting the gw.cov function to work for fit 
> points. Here's an example, using the same files...
>
>> pcavars <- names(surveypts[,4:8])
>> pcavars
> [1] "pca1" "PCA2" "PCA3" "PCA4" "PCA5"
>> surv.cov <- gw.cov(data = surveypts, vars = pcavars, fit.points = fitpix, bw = survey.bw)
> Error in gw.cov(surveypts, test, fit.points = fitpix, bw = survey.bw) :
>  unused argument(s) (fit.points = <S4 object of class "SpatialPixels">)
>
> As far as I can tell, I'm using the same syntax as in the ggwr function. 
> Thanks again. =)

?gw.cov will help. The first argument is x= not data=, and fit.points= 
should be fp=. So

surv.cov <- gw.cov(x=surveypts, vars=pcavars, fp=fitpix, bw=survey.bw)

ought to work. The functions are fairly different (no formula=), so 
assuming that the arguments have the same names suggests that actually 
reading the help pages might save time.

Please note that spgwr is only a proof of concept package and does not 
imply any affirmation of approaches using geographical weights. The 
approach may be helpful, but following Wheeler and Tiefelsdorf (2005), 
caution is recommended. Your use of PCAs may help alleviate the induced 
collinearity in GWR regression coefficients, but please always plot 
pairs() of the local coefficients!

Roger

>
> Gericke Cook
> USDA APHIS
>
>
> -----Original Message-----
>> From: modern82376 at mypacks.net
>> Sent: Mar 13, 2008 5:10 PM
>> To: r-sig-geo at stat.math.ethz.ch
>> Subject: [R-sig-Geo] writing shapefiles using write.pointShape
>>
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no