[R-sig-Geo] create ascii grid from grid.list object (fields)?

[Roger Bivand]

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