FoxHound
FoxHound

Reputation: 323

Getting NA with Akima Interpolation in R

I'm attempting to interpolate some gridded values I have to a point on a road segment using the akima packages interp function. The data points are not always all there due to the data being missing at those times. For example, see the graph below.

![Gridded Values

The lat/lng points shown are those that have dBz values from my dataset (if there is no value, our database doesn't record an entry). My goal is to use this grid to interpolate a dBz value to a specific lat/lng point (i.e. latitude = 39.95045 and longitude = -86.3574)

My data is as follows...

mrms
   tstamp: POSIXct: format: "2018-06-09 16:56:00" ...
   dBz: int 27 26 50 ...
   latitude: num 39.9 39.9 39.9 ...
   longitude: num -86.4 -86.4 -86.4

If I subset my data for a specific time (because there are multiple entries per time stamp. My data looks like the following.

        tstamp     dBz  latitude    longitude

1   6/9/2018 16:58  50  39.92501    -86.37501
2   6/9/2018 16:58  52  39.93498    -86.36501
3   6/9/2018 16:58  29  39.93498    -86.355
4   6/9/2018 16:58  38  39.91499    -86.37501
5   6/9/2018 16:58  44  39.92501    -86.36501
6   6/9/2018 16:58  23  39.92501    -86.355
7   6/9/2018 16:58  52  39.93498    -86.37501
8   6/9/2018 16:58  50  39.94498    -86.36501
9   6/9/2018 16:58  23  39.90501    -86.36501
10  6/9/2018 16:58  51  39.94498    -86.37501
11  6/9/2018 16:58  18  39.90501    -86.37501
12  6/9/2018 16:58  37  39.91499    -86.36501

After subsetting the data, I run the following command ...

interpolated <- interp(x = mrms_sub$longitude, y = mrms_sub$latitude, z = mrms_sub$dBz, xo = road$longitude, yo = road$latitude, linear = TRUE)

Which produces x = num -86.4, y = num 40, and z = num NA

I have no clue why interp keeps producing an NA value. I had this method working previously when I read all of the data in through a csv file. I've checked over datatypes for my csv method vs pulling the data in from SQL Server into a dataframe and they're the same. Any thoughts or suggestions are greatly appreciated!

UPDATE

require(tidyverse)
require(RODBC)
require(akima)

##### Road Details #####
cnx_wx <- odbcConnect("dbName", uid = "", pwd = "")
cnx_rd <- odbcConnect("dbName", uid = "", pwd = "")

queryString <- 
  "SELECT
          XDSegID,
      RoadNumber,
      RoadName,
      Bearing,
      geog.EnvelopeCenter().STAsText() as roadCenter
    FROM [itsdb1].[inrix_xd].[dbo].[__xd] as inrix_xd
    WHERE XDSegID = 1363484955 AND version = '2018-04-25'"

road <- sqlQuery(cnx_rd, queryString)

#Massage road center lat long points
road$roadCenter <- gsub("POINT|\\(|\\)", "", road$roadCenter)
road = separate(road, roadCenter, into = c("junk", "longitude", "latitude"), sep = " ")
road$longitude <- as.numeric(road$longitude)
road$latitude <- as.numeric(road$latitude)


##### Speed Data #####
queryString <- 
  "SELECT
      CAST(tstamp as smalldatetime) as tstamp,
      xdid,
      speed,
      score,
      cvalue
    FROM itsdb1.inrix_xd.dbo.xdspeeds
    WHERE tstamp >= '2018-06-09 15:00:00' AND tstamp <= '2018-06-09 18:00:00'
    AND xdid = 1363484955
    ORDER BY tstamp"

speed <- sqlQuery(cnx_rd, queryString)
#attr(speed$tstamp, "tzone") <- "GMT"
odbcClose(cnx_rd)


##### MRMS Data #####
queryString <- 
  "SELECT
      CAST(tstamp as smalldatetime) as tstamp,
      dBz,
      latitude,
      longitude
    FROM weather_mrms.dbo.v_seamless_hsr
    WHERE tstamp >= '2018-06-09 15:00:00' AND tstamp <= '2018-06-09 18:00:00'
    AND latitude BETWEEN 39.90 AND 39.95 AND longitude BETWEEN -86.38 AND -86.29
    ORDER BY tstamp"

mrms <- sqlQuery(cnx_wx, queryString)
odbcClose(cnx_wx)


##### Combined Data #####
numOfwxLevels <- nlevels(as.factor(mrms$tstamp))

i <- 1
obsTime <- mrms$tstamp[2]

output <- data.frame(0,0)
colnames(output) <- c("tstamp", "Interpolated dBz")

while (i <= numOfwxLevels) {
  mrms_sub <- filter(mrms, mrms$tstamp == obsTime)

  interpolated <- akima::interp(x = mrms_sub$longitude, y = mrms_sub$latitude, z = mrms_sub$dBz,
                         xo = road$longitude, yo = road$latitude, linear = TRUE)

  val <- c(obsTime, interpolated$z)
  newTable <- rbind(output, val)
  obsTime <- obsTime + 120
  i <- i + 1

}

Upvotes: 2

Views: 1242

Answers (1)

Anonymous coward
Anonymous coward

Reputation: 2091

Reading over the help file for it it brings up two things: 1) it recommends against colinear points, 2) it can't extrapolate using linear interpolation, which is why you're getting NA outside of your data.

If you use another interpolation it'll work, but I doubt that's what you want. You may have to look at another method.

df1terp2 <- interp(x = mrms_sub$longitude, y = mrms_sub$latitude, z = mrms_sub$dBz, xo = 39.95045, yo = -86.3574, linear = FALSE, extrap = TRUE)

$`x`
[1] 39.95045
$y
[1] -86.3574
$z
    [,1]
[1,] -78789.37

Upvotes: 2

Related Questions