On Sat, 6 Dec 2008, zack holden wrote:
Dear list, I've used R for simple statistics in the past, but this is my
first attempt at manipulating spatial data. The transition is
frustrating and I'd appreciate any advice you can give me on the problem
below.
I'm creating interpolated surfaces of PCA loadings extracted from some
precipitation data in Oregon/WA USA. I'm using the package fields. I can
create the interpolated surface, but I need to use the data from the
interpolated surface, either extracting values from the surface for
specific x,y locations, or exporting the surface to an ascii grid file,
or some format that I can either access via another R package (rgdal??)
or ArcGIS, which I've used in the past.
I'd be very grateful for suggestions and/or a snippet of example code
telling me how to make the interpolated surface created below, into a
spatial object that I can analyze further.
Thanks in advance,
Zack
Another time, please use a cluefull text-only email client - the code was
completely mixed up.
#######################################################
require(fields)
require(maptools)
require(sp)
setwd("C:/climate_data/pnw_data")
data <- readShapePoints("climate_PNW_PCLOAD.shp")
If possible use the rgdal package, it gives many more drivers for vector
and raster data - see readOGR().
# read in shapefile with X,Y and PC loadings for climate
stationsproj4string=CRS(as.character(NA)), verbose=FALSE)
No idea what this was supposed to be
# data <- as.data.frame(data)
# long <- x[,1]
# longitude column
# lat <- x[,2]
# latitude column
# minlat <- min(lat)
# maxlat <- max(lat)
# minlong <- min(long)
# maxlong <- max(long)
# grid.l<- list( X1= seq(minlong, maxlong,,100),
# X3=seq(minlat,maxlat,,100))
PredGrid <- Sobj_SpatialGrid(data)$SG
grd <- slot(PredGrid, "grid")
o <- rep(NA, prod(slot(grd, "cells.dim")))
PredGridDF <- SpatialGridDataFrame(slot(PredGrid, "grid"),
data=data.frame(o=o))
This makes an empty SpatialGridDataFrame, use the maxDim= argument to
Sobj_SpatialGrid() to control resolution.
# create grid.list object using ranges of Long. and Lat.
# make.surface.grid(grid.l)
# y <- data$pc1
# Princ. Comp. loading I want to interpolate
# LL <- cbind(long,lat)
# combine Lat and Long
# fit = Krig(LL,y, theta = 10)
fit <- Krig(coordinates(data), data$pc1, theta=10)
Use the data object directly
# fit Princ. comp. 1 to Lat long.
# surface( fit, grid.list=grid.l)
PredGridDF$pc1 <- predict(fit, coordinates(PredGridDF))
spplot(PredGridDF, "pc1", col.regions=tim.colors())
Use the predict() method on the prediction grid coordinates, assigning to
a variable in the object (treat it as you would treat a data frame).
writeGDAL(PredGridDF["pc1"], "pc1.tif", drivername="GTiff")
to write a geotiff - see other drivers with gdalDrivers().
Hope this helps,
Roger
# create interpolated surface based on Krig object
####################################################
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
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: [EMAIL PROTECTED]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
