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

[Roger Bivand]

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

> 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..

See ?tempfile, you were not using it correctly. Whether you should use it 
at all is of course also a good question - see below.

>
>> 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"))

OK so far.

>> fitpix <- SpatialPixels(fitpts, proj4string = CRS("+proj=aea +lat_0=96w 
>> +lat_1=29.5n +lat_2=45.5n"))

Why this step? It is unnecessary, just use fit.points=fitpts if the output 
Spatial*DataFrame is to be exported as a shapefile. If, on the other hand, 
you want to export it as a raster file, then don't try to write it as a 
shapefile.

>> 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"
>

If fit.points=fitpts, surv.local$SDF would be a SpatialPointsDataFrame, so 
then just say

fn <- "whatever"
writePointsShape(surv.local$SDF, fn)

and you are done. If you only want pca1, then

writePointsShape(surv.local$SDF[, "pca1"], fn)

will write just that column. If you wanted a raster, fit.points=fitpix, 
and for example:

SGDF <- as(surv.local$SDF["pca1"], "SpatialGridDataFrame")
writeAsciiGrid(SGDF, fn)

will give you an Arc ASCII grid.

>> 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"
>

You are getting exactly what you asked for - tempfile() is prepending 
"C:\temp" to the random suffix, and putting it into R's temporary 
directory as you requested.

>
> 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.

So say:

fn <- "C:/Temp/B1pca1"

and do not use tempfile(). For obvious reasons, the examples in the help 
pages use tempfile() to avoid messing up users' own work.

Roger

>
> 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
>

-- 
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