[R-sig-Geo] Question re: Intamap package / Interpolate method: Inputs NOT both Spatial Points Data Frame?

[Edzer Pebesma]

Thanks for pointing out that packages intamap and raster make
incompatible assumptions about the function (in raster: method)
interpolate. Sth we need to sort out.

In your script, you assign

gddThisYear$value <- gddThisYear[,jCtr]

it is too bad that this does not give a warning or error, apparently one
can assign S4 objects to data.frame columns -- it might even be
documented behaviour.

anyway, check the class of gddThisYear[,jCtr] yourself and try instead

gddThisYear$value <- gddThisYear[[jCtr]]

Hth,

On 10/26/2010 10:43 PM, Rick Reeves wrote:
> Hello Edzer:
> 
> (by the way, I am using R 2.12)
> 
> Yes, I had - BUT - I rechecked my script, and the 'interpolate' method
> in intamap was apparently MASKED
> by the 'interpolate' method in the subsequently loaded raster package.
> (in-progress code is below...)
> 
> I canceled the loading of the raster package, and re-ran the script with
> better results. This time,
> the line:
> 
>   interpImage <- interpolate(gddThisYear,monarchStaLocs,list(mean=TRUE,
> variance=TRUE))
> 
> results in:
> -------------------------------------------------
> R 2010-10-26 13:31:48 interpolating 133 observations, 94 prediction
> locations
> [1] "estimated time for  copula 60.5284609039258"
> Error in min(dataObs) : invalid 'type' (S4) of argument
> In addition: Warning message:
> In predictTime(nObs = dim(observations)[1], nPred =
> dim(coordinates(predictionLocations))[1],  :
> 
>  using standard model for estimating time. For better
>  platform spesific predictions, please run
>  timeModels <- generateTimeModels()
>   and save the workspace
> ------------------------------------------------
> So intamap::interpolate documentation matches the argument list. But, my
> script is still not correct.
> 
> I suspect that this error happens because I do not fully understand how
> to use 'interpolate' for this task.
> 
> My task could also be solved by an approach like that demonstrated in
> Dylan Beaudette's excellent example at:
> 
> http://casoilresource.lawr.ucdavis.edu/drupal/node/442
> 
> .. but using intamap::interpolate would be much simpler if I have made
> the correct assumptions.
> 
> What are your thoughts?
> 
> Thanks,
> Rick R
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> ##############################################################################################################
> 
> 
> interpGddData <- function()
> {
> #   require(gstat)
>    require(intamap)
> #   require(rgdal)
>    require(maptools)
> #   require(raster)
> 
>    setwd("v:/MigrationDynamics/NCEAS PROCESS")
> 
>    firstGddCol <- 5
> 
> #   setwd("/data/scientist_share/scientist/MigrationDynamics/NCEAS
> PROCESS")
> 
>    gddTable <- readShapePoints("OhioGDDStationsGP.shp")
>    monarchStaLocs <- readShapePoints("OH_Locations.shp")
> 
>    numGddCols <- 1 + (ncol(gddTable) - firstGddCol)
> 
> # The Idea: For each sampling year, generate one map for each GDD
> parameter.
> # Each map is an interpolated 'surface' from all of the GDD-calculated
> 'station'
> # locations in the input file.
> 
>    theYears <- unique(gddTable at data$YR)
>    numYears <- length(theYears)
> 
>    for (curYear in theYears)
>    {
>       gddThisYear <- gddTable[which(gddTable at data$YR==curYear),] # table
> entries: all stations, current year
> 
>    # make input table copy to append with interpolated GDD values for
> current year
> 
>       monarchStasWithGDD <- gddThisYear
> 
> # Write one output ESRI Shape File for each year, containing the Monarch
> station
> # location points with all of the interpolated GDD parameters for the
> current year appended.
> 
>       outFileName <- sprintf("OH_MonarchSta_GDD_%d.tif",curYear)
> #      theBase <- [,1:(firstGddColumn-1)]
> 
> # Create one interpolated image for each parameter, resulting in a named
> matrix file.
> 
>       for (jCtr in firstGddCol : numGddCols)
>       {
> #         fieldName <- names(gddThisYear)[jCtr] # we need this for the
> fileName
> #         fileName <- sprintf("%s_%d.tif",fieldName,curYear)
> 
>          gddThisYear$value <- gddThisYear[,jCtr]
> 
> # Note: gddThisYear needs to be upgraded to a spatialPixelsDataFrame  so
> that it is a grid
> # Despite documentation, A SpatialPointsDataFrame object is not a valid
> input.
> # first argument(spatialPoints): measurement locations (in our case, the
> GDD values)
> # second argument(spatialPixels): where to place the interpolated grid
> of measurement values.
> # I need to make an empty grid to hold the data.
> 
>          interpImage <-
> interpolate(gddThisYear,monarchStaLocs,list(mean=TRUE, variance=TRUE))
>          message("Latest interp done!")
>          browser()
> 
> # Assign the result of interpolation as an attribute to the Monarch
> Station Locations File: 'newcol'
> 
>          cbind(monarchStasWithGDD at data,newcol)
>       }
> # Write the output file
> 
>       writePointsShape(monarchStasWithGDD,outFileName)
> 
>    }
> }
> ############################################################################################################
> 
> 
> 
> 
> On 10/26/2010 1:09 PM, Edzer Pebesma wrote:
>> Are you sure you loaded package intamap before trying this?
>>
>> On 10/26/2010 09:05 PM, Rick Reeves wrote:
>>> Hello:
>>>
>>> I note a conflict in the documentation for Intamap/interpolate(),
>>> to know how to resolve:
>>>
>>> I have two SpatalPointsDataFrame objects:
>>>
>>>   growingDegreeObservations<- readShapePoints("OhioGDDStationsGP.shp")
>>> # computed surface phenomena, for individual lat/lon points
>>>   OrgStaLocs<- readShapePoints("Org_Locations.shp")                    #
>>> species observation stations, for DIFFERENT lat/lon points
>>>
>>> # Objective: Interpolate the growingDegree points into a SpatialGrid (?)
>>> raster field, and then
>>> # overlay the observation stations on to the field, and assign field
>>> values to each observation station points
>>>
>>> # according to intamap documentation, the interpolate() method accepts
>>> two SptatialPointsDataFrame objects, interpolates
>>> # the first ($value field) into a grid, and assigns interpolated
>>> observations to the predicted locations (monarchStaLocs).
>>>
>>> interpImage<- interpolate(growingDegreeObservations
>>> ,OrgStaLocs,list(mean=TRUE, variance=TRUE))
>>>
>>> However, when this line is executed, the error occurs:
>>>
>>> Error in function (classes, fdef, mtable)  :
>>>    unable to find an inherited method for function "interpolate", for
>>> signature "SpatialPointsDataFrame"
>>>
>>> I guess that this refers to the first argument, and that the argument
>>> must be a SpatialPixels or SpatialGrid data frame.
>>> Questions:
>>>
>>> 1) Do I understand the inputs and outputs of interpolate() correctly?
>>>
>>> 2) How best to transform the incorrect (SpatialPointsDataFrame) argument
>>>       to match the required inputs? E.G., create the interpolated
>>> SpatialPixelsDataFrame
>>>       of observations before calling interpolate())
>>>
>>> 3) Is there a better R technique (e.g., direct call to krige()) for
>>> computing my solution?
>>>
>>> Thanks, Rick R
>>>
>>>
>>>
>>>
>>>
>>>
>>>
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de